Práctica 1 Machine Learning 2023/2024
Master en Bioinformática, Universidad de Murcia
Fernando Lucas Ruiz (fernando.lucas@um.es)
2024-03-08
Tarea a realizar
En esta práctica se pretende analizar 3 datasets (BostonHousingData, LetterRecognition, Diabetes en el BRFSS). En cada uno de ellos vamos a responder una serie de preguntas. Se pone en práctica desde la asociación de variables pasando por PCA, el uso de algoritmos básicos como regresión lineal, logística, knn y naive bayes, y preproceso en tres datasets diferentes.
En primer lugar vamos a responder las preguntas del primer dataset BostonHousingData.
Librerias
library(tidyverse)
library(gridExtra)
library(viridis) # Para la paleta de colores viridis
library(caret)
library(kknn)
library(naivebayes)
library(pheatmap) # heatmap
library(REdaS) # Para test esfericidad de Bartlett
library(factoextra) #Gráficas PCA
library(REdaS) # Análisis PCA
library(plotly) # Gráficas PCA 3D
library(fastDummies) # codificación variables dummy
library(cvms) # matriz de confusion
library(formattable) # tablas
library(doParallel) # para paralelizar procesos
num_cores <- detectCores()
registerDoParallel(cores=num_cores)BostonHousingData
Los datos de este dataset se localizan en el paquete R mlbench.
Pregunta 1. Tamaño y tipo
¿Qué tamaño tiene? ¿De qué tipo son las variables?
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : num 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ b : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...El dataset contiene 506 filas y 14 variables.
Todas las variables son de tipo numérico excepto la variable “chas” o Charles River dummy variable.
Un primer vistazo a los datos numéricos para ver si tiene nulos, ver la distribución, buscar de variables categóricas, etc.
## BostonHousing
##
## 14 Variables 506 Observations
## --------------------------------------------------------------------------------
## crim
## n missing distinct Info Mean Gmd .05 .10
## 506 0 504 1 3.614 5.794 0.02791 0.03819
## .25 .50 .75 .90 .95
## 0.08204 0.25651 3.67708 10.75300 15.78915
##
## lowest : 0.00632 0.00906 0.01096 0.01301 0.01311
## highest: 45.7461 51.1358 67.9208 73.5341 88.9762
## --------------------------------------------------------------------------------
## zn
## n missing distinct Info Mean Gmd .05 .10
## 506 0 26 0.603 11.36 18.77 0.0 0.0
## .25 .50 .75 .90 .95
## 0.0 0.0 12.5 42.5 80.0
##
## lowest : 0 12.5 17.5 18 20 , highest: 82.5 85 90 95 100
## --------------------------------------------------------------------------------
## indus
## n missing distinct Info Mean Gmd .05 .10
## 506 0 76 0.982 11.14 7.705 2.18 2.91
## .25 .50 .75 .90 .95
## 5.19 9.69 18.10 19.58 21.89
##
## lowest : 0.46 0.74 1.21 1.22 1.25 , highest: 18.1 19.58 21.89 25.65 27.74
## --------------------------------------------------------------------------------
## chas
## n missing distinct
## 506 0 2
##
## Value 0 1
## Frequency 471 35
## Proportion 0.931 0.069
## --------------------------------------------------------------------------------
## nox
## n missing distinct Info Mean Gmd .05 .10
## 506 0 81 1 0.5547 0.1295 0.4092 0.4270
## .25 .50 .75 .90 .95
## 0.4490 0.5380 0.6240 0.7130 0.7400
##
## lowest : 0.385 0.389 0.392 0.394 0.398, highest: 0.713 0.718 0.74 0.77 0.871
## --------------------------------------------------------------------------------
## rm
## n missing distinct Info Mean Gmd .05 .10
## 506 0 446 1 6.285 0.7515 5.314 5.594
## .25 .50 .75 .90 .95
## 5.886 6.208 6.623 7.152 7.588
##
## lowest : 3.561 3.863 4.138 4.368 4.519, highest: 8.375 8.398 8.704 8.725 8.78
## --------------------------------------------------------------------------------
## age
## n missing distinct Info Mean Gmd .05 .10
## 506 0 356 0.999 68.57 31.52 17.72 26.95
## .25 .50 .75 .90 .95
## 45.02 77.50 94.07 98.80 100.00
##
## lowest : 2.9 6 6.2 6.5 6.6 , highest: 98.8 98.9 99.1 99.3 100
## --------------------------------------------------------------------------------
## dis
## n missing distinct Info Mean Gmd .05 .10
## 506 0 412 1 3.795 2.298 1.462 1.628
## .25 .50 .75 .90 .95
## 2.100 3.207 5.188 6.817 7.828
##
## lowest : 1.1296 1.137 1.1691 1.1742 1.1781
## highest: 9.2203 9.2229 10.5857 10.7103 12.1265
## --------------------------------------------------------------------------------
## rad
## n missing distinct Info Mean Gmd
## 506 0 9 0.959 9.549 8.518
##
## Value 1 2 3 4 5 6 7 8 24
## Frequency 20 24 38 110 115 26 17 24 132
## Proportion 0.040 0.047 0.075 0.217 0.227 0.051 0.034 0.047 0.261
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## tax
## n missing distinct Info Mean Gmd .05 .10
## 506 0 66 0.981 408.2 181.7 222 233
## .25 .50 .75 .90 .95
## 279 330 666 666 666
##
## lowest : 187 188 193 198 216, highest: 432 437 469 666 711
## --------------------------------------------------------------------------------
## ptratio
## n missing distinct Info Mean Gmd .05 .10
## 506 0 46 0.978 18.46 2.383 14.70 14.75
## .25 .50 .75 .90 .95
## 17.40 19.05 20.20 20.90 21.00
##
## lowest : 12.6 13 13.6 14.4 14.7, highest: 20.9 21 21.1 21.2 22
## --------------------------------------------------------------------------------
## b
## n missing distinct Info Mean Gmd .05 .10
## 506 0 357 0.986 356.7 65.5 84.59 290.27
## .25 .50 .75 .90 .95
## 375.38 391.44 396.23 396.90 396.90
##
## lowest : 0.32 2.52 2.6 3.5 3.65 , highest: 396.28 396.3 396.33 396.42 396.9
## --------------------------------------------------------------------------------
## lstat
## n missing distinct Info Mean Gmd .05 .10
## 506 0 455 1 12.65 7.881 3.708 4.680
## .25 .50 .75 .90 .95
## 6.950 11.360 16.955 23.035 26.808
##
## lowest : 1.73 1.92 1.98 2.47 2.87 , highest: 34.37 34.41 34.77 36.98 37.97
## --------------------------------------------------------------------------------
## medv
## n missing distinct Info Mean Gmd .05 .10
## 506 0 229 1 22.53 9.778 10.20 12.75
## .25 .50 .75 .90 .95
## 17.02 21.20 25.00 34.80 43.40
##
## lowest : 5 5.6 6.3 7 7.2 , highest: 46.7 48.3 48.5 48.8 50
## --------------------------------------------------------------------------------No tiene valores nulos y tienen distribuciones dispares entre variables.
Pregunta 2. Representación ejemplos
Explica qué representan los ejemplos
Este dataset incluye información sobre diferentes zonas de Boston, detalladas a través de varios parámetros. Cada registro en este dataset representa una zona de Boston, con características que influyen en la estimación del valor medio de las viviendas en esa zona. La variable dependiente es “medv”, que se refiere al valor medio de las casas habitadas por sus propietarios en miles de dólares.
Las variables del dataset son:
- crim: Tasa de criminalidad per cápita por ciudad.
- zn: Proporción de terreno residencial zonificado para lotes de más de 25,000 pies cuadrados.
- indus: Proporción de acres de negocio no minorista por ciudad.
- chas: Variable ficticia Charles River (= 1 si el tramo limita con el río; 0 de lo contrario).
- nox: Concentración de óxidos nítricos (partes por 10 millones).
- rm: Número promedio de habitaciones por vivienda.
- age: Proporción de unidades ocupadas por sus propietarios construidas antes de 1940.
- dis: Distancias ponderadas a cinco centros de empleo de Boston.
- rad: Índice de accesibilidad a autopistas radiales.
- tax: Tasa de impuesto a la propiedad de valor total por $10,000.
- ptratio: Ratio alumno-profesor por ciudad.
- b: 1000(Bk - 0.63)/^2 donde Bk es la proporción de personas negras por ciudad.
- lstat: Porcentaje de población de estatus bajo.
- medv: Valor mediano de las casas ocupadas por sus propietarios en miles de dólares.
Para ver un poco más en detalle los datos vamos a representarlos en gráficas para ver cómo son sus distribuciones o si podemos sacar alguna conclusión previa.
plots <- list()
for (i in 1:14) {
# Crea una gráfica para cada variable
p <- ggplot(BostonHousing, aes_string(x = names(BostonHousing)[i])) +
geom_density(adjust=1.5, alpha=.7, fill="#f1e6b2") +
theme_minimal() +
labs(title = names(BostonHousing)[i],
x = names(BostonHousing)[i])
plots[[colnames(BostonHousing)[i]]] <- p
}## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plots$chas <- ggplot(BostonHousing, aes(x=chas)) +
geom_bar(alpha=.7, fill="#f1e6b2", color = "black")+
theme_minimal() +
labs(title = "chas",
x = "chas")
grid.arrange(plots$crim , plots$zn, plots$indus, plots$nox, ncol=2)
Crim: Observamos que la tasa de criminalidad es baja aunque existen zonas en los datos que se llega a tasas muy altas de indice de criminalidad de hasta más de 80, que pueden ser los barrios problemáticos de la ciudad de Boston
Zn: Esta gráfica nos quiere decir que no hay muchas zonas residenciales, aunque hay dos picos intermedios sobre 23 y 80 que puede ser barrios residenciales.
Indus: Hay dos zonas de Boston con negocios de minoristas
Nox: La mayoria de zonas de Boston tienen una densidad de 0.4 y 0.6 ppm de oxido nitrico aunque hay otras menos que llegan hasta 0.9 ppm.
Rm: El número promedio de habitaciones en Boston son de aproximadamente 6 habitaciones por vivienda.
Age: Como vemos, la matoria de casas de Boston están construidas antes 1940. La cola a la izquierda nos puede sugerir que se están construyendo casas nuevas en Boston.
Dis : Esta gráfica nos indica que la mayoría de zonas se encuentran cerca de centros de empleo indicando que podría ser el centro de la ciudad, aunque la cola de la derecha nos indica que puede que los barrios residenciales se encuentren a las afueras de la ciudad.
Rad: Esta gráfica nos muestra que la mayoría de zonas se encuentran bien comunicadas con autopistas, pero hay una cierta densidad de zonas lejanas a autopistas por lo que podría indicar que son zonas residenciales o marginadas.
Tax: Esto nos indica que hay una distribución normal de las tasas excepto unas zonas que pagan más impuestos. Nos sugiere que puede ser una zona residencial de lujo.
Ptratio: La mayoría se situa sobre 20 alumnos por profesor. La distribución se ensancha en la cola izquierda porque en colegios privados el ratio suele disminuir.
B: Esta grádica nos puede sugerir que hay pocas zonas en la que la mayoría de personas son de raza caucásica.
Lstat: La tasa de personas con estatus bajo ronda el 8% pero la cola de la derecha nos sugiere que hay zonas donde el porcentaje sube por lo que puede ser zonas marginales.
Chas: La mayoría de zonas están alejadas del rio Charles
Pregunta 3. Tipo de problema
¿Es un problema de clasificación, regresión o cualquier otro? Explica por qué. Si fuera un problema de clasificación, ¿qué tienes que decir sobre la distribución de los ejemplos de cada clase? Si fuera de regresión, ¿qué tienes que decir de la distribución de valores de la variable dependiente?
En regresión, la variable dependiente o variable objetivo es continua, lo que significa que puede tomar cualquier valor dentro de un rango. En el caso de BostonHousing, la variable “medv” es un claro ejemplo de una variable continua, ya que el valor de las viviendas puede variar en un amplio espectro, teóricamente sin límites específicos dentro de los rangos observados en los datos.
El objetivo principal de un problema de regresión es predecir valores específicos de la variable dependiente basándose en una o más variables independientes (predictores). En este contexto, estamos interesados en predecir el valor de las viviendas (medv) a partir de otras características del dataset, como la tasa de criminalidad (crim), el número promedio de habitaciones por vivienda (rm), etc.
Vemos la distribución de los datos “medv” para tener una imagén inicial de esta variable.
p1 <- ggplot(BostonHousing, aes(x = medv)) +
geom_density(adjust=1.5, alpha=.7, fill="#f1e6b2") +
theme_minimal() +
labs(title = "medv", x ="") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
p2<- ggplot(BostonHousing, aes(x=medv))+
geom_boxplot(alpha=.7, fill="#f1e6b2")+
theme_minimal() +
labs(y ="")
grid.arrange(p1, p2, ncol=1)
La mayoría de casas valen alrededor de 20.000 a 23.000 USD. Vemos algunos casos interesantes que podemos tener en cuenta que hay varias casas que no siguen una tendencia normal ya que hay un pico de unas observaciones cuyos valores son de 50.000 USD.
Vamos a ver ahora cómo se correlaciona la variable “medv” con las variables predictoras
plots <- list()
for (i in 1:13) {
# Crea una gráfica para cada variable
p <- ggplot(BostonHousing, aes_string(x = names(BostonHousing)[i], y = names(BostonHousing)[14])) +
geom_point(color ="#8B1A1A", alpha = 0.7) +
labs(title = names(BostonHousing)[i],
x = names(BostonHousing)[i], y = names(BostonHousing)[14]) +
# scale_color_viridis_c(alpha = 0.7) +
xlab(names(BostonHousing)[i]) +
theme_minimal() +
ylab("medv") +
guides(color = FALSE) +
theme(axis.title.x=element_blank(),
axis.ticks.x=element_blank())
plots[[colnames(BostonHousing)[i]]] <- p
}
plots$chas <- ggplot(BostonHousing, aes(y = medv, x = chas, fill = chas)) +
geom_boxplot() +
geom_jitter(color ="#8B1A1A", size = 1, alpha = 0.7) +
theme_minimal() +
scale_fill_viridis_d() +
theme(legend.position = "none")
grid.arrange(plots$crim , plots$zn, plots$indus, plots$nox, ncol=2)
Un resumen de lo visto en los gráficos de puntos:
- Crim: el precio de los barrios con mayor indice de crimininalidad es menor.
- zn: Se aprecia una ligera relación positiva en esta variable
- indus: A mayor proporcion de negocios no minorista, el precio de la vivienda cae.
- nox: Se ve una ligera tenendecia a que cuanto mayor es la cantidad de oxido nitrico menor es el precio de las viviendas. Es decir, las zonas más alejadas de la ciudad tendrán un valor mayor en el mercado.
- rm: Se ve una clara correlación positiva entre el número de habitaciones y el precio de la vivienda.
- age: Se intuye que en los barrios más modernos el precio es ligeramente superior al casco antiguo.
- dis: Las zonas más alejadas de la ciudad parece ser que tienen los precios de viviendas más caros
- rad: Parece ser que hay barrios con poca accesibilidad a la autopista y son aquellos con un precio más bajo. Hay algunos barrios donde el precio de la vivienda es muy alto y también tiene poca accesibilidad. Puede ser que sean zonas de residencia de lujo.
- tax: Se ve una ligera correlacion negativa en esta variable.
- ptratio: Las zonas con las viviendas más caras tienen un ratio más bajo de alumno-profesorado.
- b: No se una relación clara
- lstat: Esta relación es bastante clara, ya que las personas con estatus bajo viven en zonas con el precio de la vivienda bajo
- chas: La descompesación de las clases no nos deja dar una conclusión aunque parece que las zonas cercanas al rio Charles aumentan el precio de la vivienda
correlations <- cor(BostonHousing[,c(-4)])
cor_df <- data.frame(correlations["medv", ])
cor_df$Variable <- rownames(cor_df)
cor_df <- cor_df %>% filter(Variable != "medv")
names(cor_df)[1] <- "Correlation"
cor_df <- cor_df %>% arrange(desc(Correlation))
cor_df$Variable <- factor(cor_df$Variable, levels = cor_df$Variable[order(cor_df$Correlation, decreasing = TRUE)])
ggplot(cor_df, aes(x = Variable, y = Correlation, fill= Correlation)) +
geom_col() +
coord_flip() +
labs(title = "Correlations with Medv", x = "", y = "Correlation", fill = "medv") +
scale_fill_viridis_c() 
En esta gráfica de correlación vemos que la variable mas correlacionado con la variable “medv” es aquella variable que mide el numero de personas con estatus bajo (lstat) con una correlación negativa. Luego hay otras como “ptratio”, “indus” o “tax”. La variable con mayor correlación positiva es el número de habitaciones (rm), seguida de la variable con mayor zonas residenciales (zn).
Pregunta 4. Sesgos
Supón que los posibles sesgos pudieran estar representados por los predictores del dataset. ¿Podríamos determinar visual y numéricamente, con ayuda de PCA si los hay?
El análisis de componentes principales, PCA, es una técnica estadística que permite reducir la dimensionalidad de los datos, manteniendo al mismo tiempo la mayor cantidad de variación posible. Esto se logra identificando las direcciones (componentes principales) en las cuales los datos varían más. Además de reducir la dimensionalidad, el PCA puede ofrecer insights sobre la estructura subyacente de los datos, incluidos los sesgos potenciales.
Para contrastar si nuestros datos son adecuados para hacer el análisis de componentes principales hay que pasar el test de esferecidad de Bartlett y el test de Kaiser-Meyer-Olkin.
En ambos tests rechazamos la hipotesis nula por lo que nuestros datos son adecuados para realizar el PCA Además el test KMO nos da un criterio de 0.86 que es un criterio alto de adecuación para el PCA.
## Bartlett's Test of Sphericity
##
## Call: bart_spher(x = BostonHousing[, c(-4, -14)])
##
## X2 = 4428.91
## df = 66
## p-value < 2.22e-16##
## Kaiser-Meyer-Olkin Statistics
##
## Call: KMOS(x = BostonHousing[, c(-4, -14)])
##
## Measures of Sampling Adequacy (MSA):
## crim zn indus nox rm age dis rad
## 0.9322573 0.8397090 0.8963250 0.9031321 0.7111250 0.8963153 0.8846634 0.7876067
## tax ptratio b lstat
## 0.8031895 0.8033901 0.9457321 0.8596651
##
## KMO-Criterion: 0.8560078Para aplicar el PCA de manera correcta, es esencial determinar si nuestros datos requieren un escalado apropiado.
datos_long_boston <- pivot_longer(BostonHousing, cols = -chas, names_to = "Variable", values_to = "Valor")
ggplot(datos_long_boston, aes(x = Variable, y = Valor, fill = Variable)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
theme(
legend.position = "none",
plot.title = element_text(size = 14, hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 10),
axis.title.x = element_text(margin = margin(t = 10), hjust = 0.5),
axis.title.y = element_text(margin = margin(r = 10), hjust = 0.5)
) +
ggtitle("Variables") +
xlab("Variables") +
ylab("Count")
Se debe prestar atención a variables como “tax” y “b”, las cuales se caracterizan por tener valores significativamente mayores en comparación con el resto de las variables en nuestro conjunto de datos. Este desequilibrio en la magnitud de los valores puede influir en el resultado del PCA, haciendo indispensable un ajuste de escala para asegurar una interpretación precisa y equitativa de los componentes principales.
var_exp <- pca_result_boston$sdev^2
prop_var_exp <- var_exp / sum(var_exp)
cum_var_exp <- cumsum(prop_var_exp)
df_var_exp <- data.frame(Comp = 1:length(prop_var_exp), VarExp = prop_var_exp)
df_cum_var_exp <- data.frame(Comp = 1:length(cum_var_exp), CumVarExp = cum_var_exp)
# Scree plot con los datos no escalados
ggplot(df_var_exp[1:10,], aes(x = Comp, y = VarExp)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_line(aes(group = 1), color = "blue") +
geom_point(color = "blue") +
theme_minimal() +
labs(x = "Principal components", y = "Variance", title = "Scree Plots") +
ylim(c(0,1)) +
geom_line(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), color="#8B1A1A") +
geom_point(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), color = "red") +
geom_bar(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), stat = "identity", fill = "red", alpha= 0.25) +
annotate("text", x = 9, y = 0.85, label = "Cumulative Scree \n Plot", color = "#8B1A1A", size = 4) +
geom_text(data = df_cum_var_exp[1:10,], aes(x=Comp, y = CumVarExp +0.04, label = round(CumVarExp, 2)))## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_text()`).
Las cinco componentes principales primeras explican el 85% de la varianza total, lo que nos permite avanzar con este conjunto de datos para investigar la presencia de sesgos mediante un mapa de calor que visualiza las correlaciones. Es notable que la primera componente principal por sí sola aclare más del 50% de la varianza, destacando su significativa contribución a la comprensión de la estructura de los datos.
En el siguiente heatmap vemos la correlación de las componentes principales con las variables.
sesgos <- cor(pca_result_boston$x[,1:5], BostonHousing[,c(-4,-14)])
pheatmap(sesgos,
clustering_method = "ward.D2",
color = colorRampPalette(c("#00AFBB", "white", "#FC4E07"))(255), # Paleta de colores
cellwidth = 20,
cellheight = 20,
breaks = seq(-1, 1, length.out = 256),
scale = "none",
annotation_legend = TRUE,
angle_col = 45, cutree_cols = 3)
Observamos que la componente principal 1 desempeña un papel fundamental en la diferenciación de los clusters, destacándose como el principal factor de separación. Curiosamente, las variables que presentan una menor correlación con la componente principal 1 (b, rm, zn, dis) coinciden con aquellas que muestran una mayor correlación con la variable medv, y viceversa, en el apartado de la pregunta 3 de este documento. Este patrón sugiere que dichas variables son especialmente adecuadas para inferir el precio medio de la vivienda, lo que indica su relevancia y utilidad en el análisis de los factores que determinan el valor medio de las casas en Boston.
pca <- as.data.frame(pca_result_boston$x[,1:5], stringsAsFactors=F)
pca <- cbind(BostonHousing, pca)
plots <- list()
for (i in 1:14) {
# Crea una gráfica para cada variable
p <- ggplot(pca, aes_string(x = names(pca)[15], y = names(pca)[16], color = names(pca)[i])) +
geom_point() +
labs(title = names(pca)[i],
x = names(pca)[15], y = names(pca)[16], color = names(pca)[i]) +
scale_color_gradient(low = "#D3D3D3", high = "#FC4E07") +
geom_point() +
xlab("PC 1 ") +
ylab("PC 2")+
theme_minimal()
plots[[colnames(pca)[i]]] <- p
}
plots$chas <- ggplot(pca, aes(x = PC1, y = PC2, color = chas)) +
geom_point() +
labs(title = "chas", x = "PC 1", y = "PC 2", color = "CHAS") +
scale_color_manual(values = c("0" = "#D3D3D3", "1" = "#FC4E07")) +
theme_minimal()
grid.arrange(plots$crim , plots$zn, plots$indus, plots$nox, ncol=2)
Como ya veiamos en el heatmap, hay un sesgo claro en nuestros datos ya que nuestros datos del PCA generan dos clusters grandes. Observamos que seis barrios se distancian significativamente del resto, ubicados en la zona inferior derecha del gráfico.
El cluster con forma redondeada, que podría interpretarse como el núcleo urbano de Boston, se caracteriza por su óptima accesibilidad a las autopistas (rad), un régimen tributario más elevado (tax), una concentración superior de edificaciones antiguas (age), niveles más altos de óxidos de nitrógeno (nox), y una mayor densidad de zonas industriales (indus).
Por contraste, el segundo cluster alargado sugiere una transición hacia la periferia, como lo indica el aumento en la distancia media a centros de empleo (dis). Este patrón se ve reflejado en la reducción progresiva de zonas industriales (indus) hacia las afueras, dando paso a un incremento de áreas residenciales (zn), especialmente evidente en la extensión del cluster elongado. Adicionalmente, la presencia de viviendas con un mayor número promedio de habitaciones (rm) hacia la izquierda del cluster refuerza esta noción de expansión urbana.
Un detalle revelador es la existencia de un barrio céntrico con un porcentaje inusualmente bajo de residentes negros (b). Sin embargo, los datos no proporcionan más características distintivas de este barrio.
La proximidad al río Charles (chas) en el cluster redondeado parece delinear una zona de exclusividad en el corazón de la ciudad, donde los valores inmobiliarios (medv) son sustancialmente más altos y las viviendas ostentan un número mayor de habitaciones en comparación con el resto del centro urbano.
Enfocando la atención en los seis outliers en la esquina inferior derecha, se aprecian índices de criminalidad (crim) más elevados y un porcentaje más alto de personas de estatus socioeconómico bajo (lstat). Estos barrios presentan una marcada diversidad racial ya que, o son exlusivamente personas negras o personas blancas (b), y están situados cerca de las autopistas, caracteristica de barrios marginales por el tráfico de drogas.
Finalmente, al examinar la variable dependiente medv, se deduce que la distancia del núcleo urbano correlaciona con un incremento en el valor de las viviendas, que tienden a ser más espaciosas y situadas en zonas predominantemente residenciales. El sector adyacente al río Charles en el núcleo urbano se destaca por su alto valor inmobiliario. En contraste, las áreas con una tasa de criminalidad más alta tienden a coincidir con los precios de vivienda más bajos, lo que señala una relación inversa entre seguridad y valor inmobiliario.
Vamos a ver ahora la importancia de las variables a cada una de las componentes principales.
p1 <- fviz_cos2(pca_result_boston, choice = "var", axes = 1, top= 5, fill = "#00AFBB", ggtheme = theme_minimal()) + labs(y = "Cos2", title = "Dim 1")
p2 <- fviz_cos2(pca_result_boston, choice = "var", axes = 2, top= 5, fill = "#00AFBB", ggtheme = theme_minimal()) + labs(y = "Cos2", title = "Dim 2")
p3 <- fviz_cos2(pca_result_boston, choice = "var", axes = 3, top= 5, fill = "#00AFBB", ggtheme = theme_minimal()) + labs(y = "Cos2", title = "Dim 3")
p4 <- fviz_cos2(pca_result_boston, choice = "var", axes = 4, top= 5, fill = "#00AFBB", ggtheme = theme_minimal()) + labs(y = "Cos2", title = "Dim 4")
grid.arrange(p1, p2, p3, p4, ncol=2)
El cuadrado del coseno (Cos2) de una variable en una componente principal indica cuánto contribuye esa variable a la varianza explicada por esa componente principal. Para cada componente principal, se presentan las cinco variables con el mayor Cos2 en esa dimensión particular.
Observamos que la primera y segunda componentes principales tiene varias variables que explican su varianza. Aunque la contribución de las variables a la varianza de la componente principal 2 es notable, es menos significativa que las contribuciones de las variables a la componente principal 1, como se puede ver por los valores más bajos de Cos2.
En la componente principal 3, la variable rm tiene un Cos2 mucho más alto que las demás, indicando que es un factor más dominante en esta componente.
En la componente 4 ptratio tiene la contribución más alta, pero todas tienen valores relativamente bajos, lo que sugiere que la componente principal 4 captura menos varianza en las variables originales que las componentes anteriores.
pca_rotation <- round(data.frame(pca_result_boston$rotation[,1:4]),2)
pca_rotation <- pca_rotation[order(-pca_rotation[,1]),]
formattable(pca_rotation,
list(
area(col= PC1:PC4) ~ formatter("span", style = x ~ ifelse(x > 0,
style(color = "green", font.weight = "bold"), style(color = "red", font.weight = "bold")))
))| PC1 | PC2 | PC3 | PC4 | |
|---|---|---|---|---|
| indus | 0.35 | 0.11 | 0.03 | -0.01 |
| nox | 0.34 | 0.17 | 0.24 | 0.15 |
| tax | 0.34 | -0.32 | 0.08 | -0.14 |
| rad | 0.32 | -0.38 | 0.11 | -0.22 |
| age | 0.31 | 0.32 | 0.16 | 0.03 |
| lstat | 0.31 | 0.03 | -0.30 | 0.39 |
| crim | 0.25 | -0.40 | 0.07 | 0.08 |
| ptratio | 0.21 | -0.17 | -0.49 | -0.61 |
| rm | -0.19 | -0.08 | 0.68 | -0.39 |
| b | -0.20 | 0.34 | -0.19 | -0.36 |
| zn | -0.26 | -0.44 | 0.09 | 0.30 |
| dis | -0.32 | -0.33 | -0.25 | 0.08 |
Los vectores de carga de un PCA representan la contribución de cada variable original a las componentes principales que se han calculado. Cada columna en los datos que has proporcionado representa una componente principal, y cada fila representa una variable. El valor de cada celda en la matriz de carga es el coeficiente de la variable original en la componente principal correspondiente. Estos coeficientes pueden ser interpretados como la importancia o peso de esa variable para esa componente principal específica. Un valor alto (positivo o negativo) indica que la variable tiene una fuerte relación con esa componente principal, mientras que un valor cercano a cero indica que la variable tiene poco o ningún efecto en esa dimensión. Los signos de los valores de carga indican la dirección de la relación entre la variable y la componente principal, pero en el contexto del PCA, solo la magnitud de estos valores es relevante para determinar la importancia de la variable. Los signos son arbitrarios y pueden cambiar si se ejecuta el PCA nuevamente debido a la inversión del eje.
En la componente principal 1, las tendencias generales observadas en esta componente principal son valores positivos altos para “indus”, “nox”, “age”, “rad”, y “tax”, lo cual sugiere que estas variables contribuyen significativamente y en la misma dirección a la primera componente principal. Estas variables podrían estar asociadas con el centro urbano o el desarrollo industrial.
En la componente principal 2, las variables “zn” y “crim” tienen los valores de carga más altos por magnitud pero en direcciones opuestas. “zn” representa las zonas residenciales, mientras que “crim” podría ser más alta en áreas urbanas densamente pobladas.
La tercera componente principal está dominada por “rm” con un valor alto y positivo, lo que significa que “rm” es el contribuyente principal a esta componente.
Para la cuarta componente principal, “ptratio” y “b” tienen los valores más altos, pero ptratio tiene una carga negativa fuerte, lo que podría indicar una relación inversa con esta componente. Esta dimensión podría capturar aspectos del entorno social o educativo en las áreas del estudio.
A continuación, muestro gráficamente los resultados de las dos primeras componentes principapeles para ver la dirección e importancia de las mismas con lo descrito en este apartado.
fviz_pca_biplot(pca_result_boston,
geom.ind = "point",
geom_var = c("arrow", "text"),
col.var = "cos2",
col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = F,
)
Pregunta 5. Elección algoritmo
¿Qué algoritmo de entre los que conoces aplicarías según el problema para el que lo has identificado?
Dada la naturaleza de los datos de Bostonhousing, considero que se debe utilizar un algoritmo de regresión para predecir la variable medv con respecto a las demás variables. Debido al tamaño de los datos los mejores algoritmos para esta práctica es algoritmo de k-vecinos más cercanos (KNN).
KNN es un algoritmo no paramétrico y puede capturar relaciones no lineales entre las variables sin una formulación específica de un modelo. KNN puede ser útil para identificar patrones locales y no lineales que la regresión lineal podría pasar por alto, especialmente en áreas donde las propiedades tienen características distintivas que solo se aplican a un vecindario específico.
El algoritmo KNN es intuitivo y fácil de entender ya que predice el valor de una nueva observación promediando los valores de los K vecinos más cercanos. KNN tiene en cuenta la localidad de los datos, lo que podría ser particularmente relevante para datos inmobiliarios donde la ubicación y el vecindario juegan un papel crucial en el precio de las viviendas.
Pregunta 6. Transformación
¿Sería necesario aplicar alguna transformación a los datos? ¿Por qué? Si la respuesta es afirmativa, explica qué tendría de beneficioso hacerlo.
Es importante prestar atención al preprocesamiento de los datos al utilizar KNN. Dado que KNN depende de la distancia entre las observaciones, las variables deben ser normalizadas o estandarizadas para que tengan la misma escala. De lo contrario, las variables con rangos más grandes podrían influir desproporcionadamente en la predicción.
También debemos pasar la variable categórica chas a numérica mediante codificación Dummy. En este caso, como solamente hay dos niveles, sería solamente pasar la variable a numerica entre 0 y 1 sin añadir columnas adicionales.
BostonHousing.scaled <- data.frame(scale(BostonHousing[,-4]))
BostonHousing.scaled <- cbind(BostonHousing.scaled, chas = BostonHousing$chas)
BostonHousing.scaled$chas <- as.numeric(BostonHousing.scaled$chas) - 1
datos_long_boston <- pivot_longer(BostonHousing.scaled, cols = -chas, names_to = "Variable", values_to = "Valor")
ggplot(datos_long_boston, aes(x = Variable, y = Valor, fill = Variable)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
theme(
legend.position = "none",
plot.title = element_text(size = 14, hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 10),
axis.title.x = element_text(margin = margin(t = 10), hjust = 0.5),
axis.title.y = element_text(margin = margin(r = 10), hjust = 0.5)
) +
ggtitle("Variables") +
xlab("Variables") +
ylab("Count")
Pregunta 7. Algoritmo
Aplica el algoritmo designado en la pregunta 5 y optimiza el estadístico o estadísticos que consideres pertinente para medir su rendimiento mediante la selección de un modelo óptimo. Comenta los resultados
En primer, lugar vamos a generar los datos de entrenamiento y test de los datos de BostonHousing escalados para entrenar el algoritmo de KNN.
set.seed(1234)
boston.TrainIdx.80<- createDataPartition(BostonHousing.scaled$medv,
p=0.8, #Genera un 80% para train, 20% para test
list = FALSE, # resultados en una matriz
times = 1) #Genera solamente una partición 80/20
trainSet.boston <- BostonHousing.scaled[boston.TrainIdx.80, ]
testSet.boston <- BostonHousing.scaled[-boston.TrainIdx.80, ]La función kknn es parte del paquete kknn, que está integrado con caret. Algunos de los hiperparámetros que se pueden ajustar al utilizar la función kknn a través de caret incluyen:
kmax: El número de vecinos más cercanos a considerar en la votación o el promedio. Es el hiperparámetro principal del KNN. Un k pequeño puede hacer que el modelo sea sensible al ruido de los datos (sobreajuste), mientras que un k grande puede hacer que el modelo sea demasiado general (subajuste).
distance: Define la métrica de distancia a utilizar para calcular la cercanía entre puntos. Por ejemplo, una distancia de 2 se refiere a la distancia euclidiana, mientras que 1 se refiere a la distancia de Manhattan. Distancias fraccionarias también son posibles y pueden ser útiles para capturar estructuras no lineales.
kernel: Especifica la función de ponderación a aplicar a los vecinos. Algunas opciones incluyen “rectangular” (todos los vecinos ponderan igual), “triangular”, “epanechnikov”, “biweight”, “triweight”, “cos”, entre otros. La ponderación puede afectar cómo influyen los vecinos en la predicción final, dando más peso a los vecinos más cercanos.
Es importante destacar que el ajuste de hiperparámetros debe hacerse cuidadosamente. Por ejemplo, seleccionar k puede requerir la implementación de técnicas como la validación cruzada para encontrar un valor que equilibre bien el sesgo y la varianza. La elección de la métrica de distancia y la función de kernel debe basarse en el conocimiento del problema y el dominio de aplicación, así como en la experimentación y validación.
Debido a que el numero de muestras que tenemos no es muy grande, he decidido optar por un entrenamiento de validación cruzada repetitiva con 3 repeticiones y 10 particiones. Permito el análisis en paralelo para aumentar el rendimiento
set.seed(1234)
kknnControl.boston <- trainControl(method = "repeatedcv",
number = 10,
repeats = 3,
seeds = NULL,
returnResamp = "final",
allowParallel = T)En un primer intento de buscar el mejor modelo realizo una parrilla de hiperparámetros con los siguientes valores dando 21 posibilidades combinatorias:
kmax: 5, 7, 9 , 11
distance: 1, 2, 3
kernel: rectangular, triangular, optimal
set.seed(1234)
y = trainSet.boston$medv
x = trainSet.boston[,-13]
mygrid.boston = expand.grid(kmax = seq(from=5,to=11,2),
distance = seq(1, 3, 1),
kernel = c("rectangular","triangular","optimal"
# ,"epanechnikov", "biweight", "triweight"
# ,"cos", "inv", "gaussian","rank"
))
myfit.boston.cv10.rp3 <- train(y=y, x = x,
method = "kknn",
metric = "Rsquared",
trControl = kknnControl.boston,
tuneGrid=mygrid.boston
)
(myfit.boston.cv10.rp3.best <- subset(myfit.boston.cv10.rp3$results, kmax == myfit.boston.cv10.rp3$bestTune$kmax & distance == myfit.boston.cv10.rp3$bestTune$distance & kernel == myfit.boston.cv10.rp3$bestTune$kernel))
En este primer modelo vemos que los mejores hiperparámetros son con una kmax de 5, distancia Manhattan de 1 y el kernel optimal. La RMSE y Rsquared dan valores muy buenos para este modelo.
Vamos a probar con otros tipos de kernel para ver si mejoramos el modelo. Los hiperparámetros son:
kmax: 5, 7, 9 , 11
distance: 1
kernel: optimal, epanechnikov, biweight, triweight, cos, gaussian, rank
set.seed(1234)
mygrid.boston = expand.grid(kmax = seq(from=5,to=11,3),
distance = 1,
kernel = c(#"rectangular","triangular",
"optimal"
,"epanechnikov", "biweight", "triweight"
,"cos", "gaussian","rank"
))
myfit.boston.cv10.rp3.2 <- train(y=y, x = x,
method = "kknn",
metric = "Rsquared",
trControl = kknnControl.boston,
tuneGrid=mygrid.boston,
#preProcess = c("center", "scale")
)
(myfit.boston.cv10.rp3.2best <- subset(myfit.boston.cv10.rp3.2$results, kmax == myfit.boston.cv10.rp3.2$bestTune$kmax & distance == myfit.boston.cv10.rp3.2$bestTune$distance & kernel == myfit.boston.cv10.rp3.2$bestTune$kernel))
El mejor modelo posible lo seguimos encontrando con los hiperparámetros anteriores, por lo que no mejoramos el modelo. Como vemos en las siguiente gráficas, casi todas las iteraciones de nuestro modelo de entrenamiento nos da como resultado unas Rsquared de entorno 0.85 y un RMSE de 0.39 sin grandes cambios, aunque hay unas pocas iteraciones con un RMSE un poco más alto y un Rsquared más bajo.
resample <- myfit.boston.cv10.rp3.2$resample
p1 <- ggplot(data = resample, aes(x = RMSE)) +
geom_density(size = 2,color = "#8B1A1A") +
theme_minimal() +
labs(title = "RMSE",
x = "RMSE")+
geom_vline(xintercept = median(resample$RMSE),
linetype = "dashed") +
annotate("text", x = 0.5, y = 4,
label = paste("RMSE", "\n", round(median(resample$RMSE), 2)), color = "#8B1A1A", size = 4)
p2 <- ggplot(data = resample, aes(x = Rsquared)) +
geom_density(size = 2, color = "#f1e6b2") +
theme_minimal() +
labs(title = "Rsquared",
x = "Rsquared")+
geom_vline(xintercept = median(resample$Rsquared),
linetype = "dashed") +
annotate("text", x = 0.75, y = 6,
label = paste("Rsquared", "\n", round(median(resample$Rsquared), 2)), color = "#c0b88e", size = 4)
grid.arrange(p1, p2, ncol = 2)
En la siguiente tabla vemos las variables más importantes para nuestro modelo situando a “lstat” y “rm” como las variables más importantes. “Chas” no influye nada en nuestro modelo ya que los datos están muy desbalanceados.
var.import.boston <- varImp(myfit.boston.cv10.rp3.2)
var.import.boston.ordered <- var.import.boston$importance %>%
arrange(desc(Overall))
formattable(var.import.boston.ordered,
list(
Overall = color_tile("white", "#FC4E07"))
)| Overall | |
|---|---|
| lstat | 100.00000 |
| rm | 85.63170 |
| ptratio | 34.51859 |
| indus | 30.43361 |
| crim | 28.26487 |
| tax | 27.10310 |
| nox | 23.96533 |
| b | 19.83734 |
| age | 18.52760 |
| rad | 17.22180 |
| zn | 16.61144 |
| dis | 12.80692 |
| chas | 0.00000 |
Vemos que los valores predecidos caen en lugares parecidos a los valores de train y de test sin que haya ningún outlier aparente.
preds.boston <- predict(myfit.boston.cv10.rp3.2$finalModel, newdata = testSet.boston)
df.train <- trainSet.boston
df.test <- testSet.boston
df.train$trainORtest <- "train"
df.test$trainORtest <- "test"
boston_df_plot.boston <- rbind(df.train, df.test)
prediction_test.boston <- cbind(testSet.boston, preds.boston)
plots <- list()
for (i in 1:13){
col_name <- names(boston_df_plot.boston)[i]
p <- ggplot(boston_df_plot.boston, aes_string(x = col_name, y = "medv")) +
geom_point(aes(color = trainORtest)) +
geom_point(data = prediction_test.boston, aes_string(x = col_name, y = "preds.boston", color = "'Predicciones'"), show.legend = TRUE) +
scale_color_manual(values = c("train" = "#f1e6b2", "test" = "#b2d8b2", "Predicciones" = "black"),
name = "", labels = c("Predicciones", "Test", "Train")) +
theme_minimal()
plots[[col_name]] <- p
}
plots$chas <- ggplot(boston_df_plot.boston, aes(y = medv, x = chas, color = trainORtest)) +
geom_jitter() +
geom_jitter(data = prediction_test.boston, aes(x = chas, y = preds.boston, color = "Predictions"), show.legend = TRUE) +
scale_color_manual(values = c("train" = "#f1e6b2", "test" = "#b2d8b2", "Predictions" = "black"),
name = "", labels = c("Predictions", "Test", "Train")) +
theme_minimal() +
geom_vline(xintercept = 0.5, color = "#8B1A1A", linetype = "dashed")
grid.arrange(plots$crim , plots$zn, plots$indus, plots$nox, ncol=2)
El análisis del error de predicción en relación con los valores reales de la variable “medv” en el conjunto de test revela una tendencia diferenciada en función del rango de precios de las viviendas. Observamos que para las viviendas con un valor cercano a la mediana de medv, las predicciones son notablemente precisas, reflejando un error bajo. No obstante, para las viviendas con un valor de “medv” superior al promedio, el modelo tiende a tener un rendimiento menos exacto, resultando en errores de predicción más elevados.
testSet.boston %>%
cbind(preds.boston) %>%
mutate(RMSE = sqrt((preds.boston - medv)^2)) %>%
ggplot(aes(x = medv, y = RMSE)) +
geom_line(color ="#8B1A1A") +
geom_point(color ="#c0b88e") +
theme_minimal() +
labs(title = "Predictive error test data", x = "Medv test data")
Finalmente, para ver si incurrimos en overfitting o underfitting, calculamos el RMSE de los valores predichos con los datos de entrenamiento y con lo datos de test. Observamos que en los datos de entrenamiento el error es menor que en los datos de test pero la diferencia no es muy grande, siendo el error muy bajo.
Encontramos un modelo bueno para predecir el precio de la vivienda con estas variables en la ciudad de Boston.
preds.boston.train <- predict(myfit.boston.cv10.rp3.2$finalModel, newdata = trainSet.boston)
cat("El error RMSE de los datos de entrenamiento es",sqrt((1/nrow(trainSet.boston)) * sum((trainSet.boston$medv - preds.boston.train)^2)),"\n")## El error RMSE de los datos de entrenamiento es 0.2160162cat("El error RMSE de los datos de evaluación es", sqrt((1/nrow(testSet.boston)) * sum((testSet.boston$medv - preds.boston)^2)),"\n")## El error RMSE de los datos de evaluación es 0.374676LetterRecognition
Los datos de este dataset se localizan en el paquete R mlbench.
Pregunta 1. Tamaño y tipo
¿Qué tamaño tiene? ¿De qué tipo son las variables?
## 'data.frame': 20000 obs. of 17 variables:
## $ lettr: Factor w/ 26 levels "A","B","C","D",..: 20 9 4 14 7 19 2 1 10 13 ...
## $ x.box: num 2 5 4 7 2 4 4 1 2 11 ...
## $ y.box: num 8 12 11 11 1 11 2 1 2 15 ...
## $ width: num 3 3 6 6 3 5 5 3 4 13 ...
## $ high : num 5 7 8 6 1 8 4 2 4 9 ...
## $ onpix: num 1 2 6 3 1 3 4 1 2 7 ...
## $ x.bar: num 8 10 10 5 8 8 8 8 10 13 ...
## $ y.bar: num 13 5 6 9 6 8 7 2 6 2 ...
## $ x2bar: num 0 5 2 4 6 6 6 2 2 6 ...
## $ y2bar: num 6 4 6 6 6 9 6 2 6 2 ...
## $ xybar: num 6 13 10 4 6 5 7 8 12 12 ...
## $ x2ybr: num 10 3 3 4 5 6 6 2 4 1 ...
## $ xy2br: num 8 9 7 10 9 6 6 8 8 9 ...
## $ x.ege: num 0 2 3 6 1 0 2 1 1 8 ...
## $ xegvy: num 8 8 7 10 7 8 8 6 6 1 ...
## $ y.ege: num 0 4 3 2 5 9 7 2 1 1 ...
## $ yegvx: num 8 10 9 8 10 7 10 7 7 8 ...El dataset contiene 20.000 filas y 17 variables.
Todas las variables son de tipo numérico excepto la variable lettr que es la variable dependiente.
Un primer vistazo a los datos numéricos para ver si tiene nulos, la distribución, búsqueda de variables categóricas, etc
## LetterRecognition
##
## 17 Variables 20000 Observations
## --------------------------------------------------------------------------------
## lettr
## n missing distinct
## 20000 0 26
##
## lowest : A B C D E, highest: V W X Y Z
## --------------------------------------------------------------------------------
## x.box
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.972 4.024 2.093 1 2
## .25 .50 .75 .90 .95
## 3 4 5 7 7
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 132 1261 2909 4157 4477 3169 1894 1006 513 284 121
## Proportion 0.007 0.063 0.145 0.208 0.224 0.158 0.095 0.050 0.026 0.014 0.006
##
## Value 11 12 13 14 15
## Frequency 48 20 4 3 2
## Proportion 0.002 0.001 0.000 0.000 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## y.box
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.991 7.035 3.73 1 2
## .25 .50 .75 .90 .95
## 5 7 9 11 12
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 709 778 530 1330 1342 1566 1705 2302 2180 2702 2211
## Proportion 0.035 0.039 0.026 0.066 0.067 0.078 0.085 0.115 0.109 0.135 0.111
##
## Value 11 12 13 14 15
## Frequency 1625 321 271 198 230
## Proportion 0.081 0.016 0.014 0.010 0.012
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## width
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.975 5.122 2.23 2 3
## .25 .50 .75 .90 .95
## 4 5 6 8 9
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 195 385 1285 1994 3816 4262 3641 1946 1418 679 237
## Proportion 0.010 0.019 0.064 0.100 0.191 0.213 0.182 0.097 0.071 0.034 0.012
##
## Value 11 12 13 14 15
## Frequency 91 39 6 4 2
## Proportion 0.005 0.002 0.000 0.000 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## high
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.98 5.372 2.535 1 2
## .25 .50 .75 .90 .95
## 4 6 7 8 8
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 365 883 1304 1559 2718 2675 3656 2695 3559 347 103
## Proportion 0.018 0.044 0.065 0.078 0.136 0.134 0.183 0.135 0.178 0.017 0.005
##
## Value 11 12 13 14 15
## Frequency 76 31 10 15 4
## Proportion 0.004 0.002 0.000 0.001 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## onpix
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.976 3.506 2.377 1 1
## .25 .50 .75 .90 .95
## 2 3 5 6 8
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 830 2437 4153 3939 3157 2153 1379 857 519 283 142
## Proportion 0.042 0.122 0.208 0.197 0.158 0.108 0.069 0.043 0.026 0.014 0.007
##
## Value 11 12 13 14 15
## Frequency 85 40 13 7 6
## Proportion 0.004 0.002 0.001 0.000 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## x.bar
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.96 6.898 2.156 3 4
## .25 .50 .75 .90 .95
## 6 7 8 9 10
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 148 179 167 680 1069 1780 2717 6024 4019 1752 802
## Proportion 0.007 0.009 0.008 0.034 0.053 0.089 0.136 0.301 0.201 0.088 0.040
##
## Value 11 12 13 14 15
## Frequency 338 201 67 43 14
## Proportion 0.017 0.010 0.003 0.002 0.001
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## y.bar
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.963 7.5 2.488 3 5
## .25 .50 .75 .90 .95
## 6 7 9 11 12
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 47 134 393 488 714 877 2506 6010 3753 1599 1252
## Proportion 0.002 0.007 0.020 0.024 0.036 0.044 0.125 0.300 0.188 0.080 0.063
##
## Value 11 12 13 14 15
## Frequency 1131 663 192 165 76
## Proportion 0.057 0.033 0.010 0.008 0.004
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## x2bar
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.983 4.629 2.878 1 2
## .25 .50 .75 .90 .95
## 3 4 6 8 9
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 422 1084 2693 3424 3199 2982 2340 1422 1013 436 235
## Proportion 0.021 0.054 0.135 0.171 0.160 0.149 0.117 0.071 0.051 0.022 0.012
##
## Value 11 12 13 14 15
## Frequency 150 158 120 184 138
## Proportion 0.007 0.008 0.006 0.009 0.007
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## y2bar
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.984 5.179 2.66 1 2
## .25 .50 .75 .90 .95
## 4 5 7 8 9
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 269 822 1852 1914 3011 3243 3099 2623 1688 874 318
## Proportion 0.013 0.041 0.093 0.096 0.151 0.162 0.155 0.131 0.084 0.044 0.016
##
## Value 11 12 13 14 15
## Frequency 128 48 31 46 34
## Proportion 0.006 0.002 0.002 0.002 0.002
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## xybar
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.971 8.282 2.736 5 6
## .25 .50 .75 .90 .95
## 7 8 10 12 13
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 145 62 97 90 267 819 2751 5578 1668 1971 2563
## Proportion 0.007 0.003 0.005 0.004 0.013 0.041 0.138 0.279 0.083 0.099 0.128
##
## Value 11 12 13 14 15
## Frequency 1799 1043 772 279 96
## Proportion 0.090 0.052 0.039 0.014 0.005
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## x2ybr
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.972 6.454 2.876 2 3
## .25 .50 .75 .90 .95
## 5 6 8 10 11
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 188 430 887 726 1495 2457 5666 2598 1495 1408 888
## Proportion 0.009 0.021 0.044 0.036 0.075 0.123 0.283 0.130 0.075 0.070 0.044
##
## Value 11 12 13 14 15
## Frequency 935 472 200 135 20
## Proportion 0.047 0.024 0.010 0.007 0.001
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## xy2br
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.957 7.929 2.232 4 5
## .25 .50 .75 .90 .95
## 7 8 9 10 12
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 1 9 67 269 701 1335 1874 2807 6548 3000 1437
## Proportion 0.000 0.000 0.003 0.013 0.035 0.067 0.094 0.140 0.327 0.150 0.072
##
## Value 11 12 13 14 15
## Frequency 926 409 313 219 85
## Proportion 0.046 0.020 0.016 0.011 0.004
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## x.ege
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.972 3.046 2.502 0 0
## .25 .50 .75 .90 .95
## 1 3 4 6 8
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 2461 2568 4213 4779 1500 1363 1264 722 587 246 154
## Proportion 0.123 0.128 0.211 0.239 0.075 0.068 0.063 0.036 0.029 0.012 0.008
##
## Value 11 12 13 14 15
## Frequency 81 29 17 12 4
## Proportion 0.004 0.001 0.001 0.001 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## xegvy
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.935 8.339 1.664 6 6
## .25 .50 .75 .90 .95
## 8 8 9 10 11
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 1 13 17 34 79 416 1592 2516 7624 3437 2394
## Proportion 0.000 0.001 0.001 0.002 0.004 0.021 0.080 0.126 0.381 0.172 0.120
##
## Value 11 12 13 14 15
## Frequency 1437 348 72 16 4
## Proportion 0.072 0.017 0.004 0.001 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## y.ege
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.986 3.692 2.89 0 0
## .25 .50 .75 .90 .95
## 2 3 5 7 8
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 2472 2040 2475 3078 3091 2048 1723 1227 973 613 154
## Proportion 0.124 0.102 0.124 0.154 0.155 0.102 0.086 0.061 0.049 0.031 0.008
##
## Value 11 12 13 14 15
## Frequency 64 22 13 3 4
## Proportion 0.003 0.001 0.001 0.000 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## yegvx
## n missing distinct Info Mean Gmd .05 .10
## 20000 0 16 0.927 7.801 1.714 5 6
## .25 .50 .75 .90 .95
## 7 8 9 10 11
##
## Value 0 1 2 3 4 5 6 7 8 9 10
## Frequency 2 17 30 130 478 992 1827 3472 8047 2358 1578
## Proportion 0.000 0.001 0.002 0.006 0.024 0.050 0.091 0.174 0.402 0.118 0.079
##
## Value 11 12 13 14 15
## Frequency 868 137 49 13 2
## Proportion 0.043 0.007 0.002 0.001 0.000
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------No tiene valores nulos y las distribuciones van de 0 a 15 en todas las variables menos en la variable dependiente lettr.
Pregunta 2. Representación ejemplos
Explica qué representan los ejemplos
Los datos de LetterRecognition del paquete mlbench en R están diseñados para la tarea de identificar letras del alfabeto inglés en mayúsculas a partir de imágenes en píxeles blanco y negro. Cada imagen de letra se basa en 20 fuentes diferentes, y cada letra de estas fuentes se distorsionó aleatoriamente para producir un conjunto de 20,000 estímulos únicos. Estos estímulos se convirtieron en 16 atributos numéricos primitivos, como momentos estadísticos y recuentos de bordes, que se escalaron para ajustarse a un rango de valores enteros de 0 a 15.
El dataset contiene 20,000 observaciones con 17 variables. La primera variable es categórica y tiene niveles de la A a la Z, representando cada una de las letras del alfabeto. Las 16 variables restantes son numéricas y representan diversas características estadísticas y geométricas de las imágenes de las letras.
x.box: La posición horizontal de la caja que encierra la letra (el rectángulo más pequeño que puede contener la imagen de la letra).
y.box: La posición vertical de la caja que encierra la letra.
width: El ancho de la caja en píxeles.
high: La altura de la caja en píxeles.
onpix: El total de píxeles encendidos (es decir, píxeles que contienen parte de la letra) en la imagen.
x.bar: La media de la posición horizontal de todos los píxeles encendidos.
y.bar: La media de la posición vertical de todos los píxeles encendidos.
x2bar: La varianza de la posición horizontal de todos los píxeles encendidos.
y2bar: La varianza de la posición vertical de todos los píxeles encendidos.
xybar: La correlación entre la posición horizontal y vertical de los píxeles encendidos.
x2ybr: La media del producto de x por la posición vertical de los píxeles encendidos.
xy2br: La media del producto de y por la posición horizontal de los píxeles encendidos.
x.ege: La cantidad de bordes horizontales (transiciones de claro a oscuro).
xegvy: La correlación de x.ege con la variable y.
y.ege: La cantidad de bordes verticales (transiciones de claro a oscuro).
yegvx: La correlación de y.ege con la variable x.
Pregunta 3. Tipo de problema
¿Es un problema de clasificación, regresión o cualquier otro? Explica por qué. Si fuera un problema de clasificación, ¿qué tienes que decir sobre la distribución de los ejemplos de cada clase? Si fuera de regresión, ¿qué tienes que decir de la distribución de valores de la variable dependiente?
Este conjunto de datos es un problema de clasificación. La razón es que el objetivo es categorizar cada representación digital de una letra en una de varias clases predefinidas, que son las letras del alfabeto. Por lo tanto, el conjunto de datos es un ejemplo clásico de un problema de clasificación multiclase, donde el modelo de aprendizaje automático debe predecir cuál de las múltiples categorías (lettr) corresponde a las características dadas de una imagen.
Para un problema de clasificación, es crucial examinar la distribución de los ejemplos entre las clases. Una distribución equilibrada de ejemplos asegura que el modelo de aprendizaje automático tenga suficientes datos para aprender las características distintivas de cada clase. Por otro lado, una distribución desequilibrada puede llevar a un modelo sesgado que prefiera las clases con más ejemplos.
A continuación mostramos las distribuciones de las frecuencias de cada clase de la variable “lettr”. Vemos que la distribución es bastante homogénea entre clases ya que ronda las 750 instancias de cada letra.
##
## A B C D E F G H I J K L M N O P Q R S T
## 789 766 736 805 768 775 773 734 755 747 739 761 792 783 753 803 783 758 748 796
## U V W X Y Z
## 813 764 752 787 786 734LetterRecognition %>%
group_by(lettr) %>%
summarize(Frequency = n()) %>%
ungroup() %>%
ggplot(aes(x = lettr, y = Frequency)) +
geom_segment(aes(x = lettr, xend = lettr, y = 0, yend = Frequency)) +
geom_point(size = 5, aes(fill = lettr), shape = 21, stroke = 2, alpha = 0.7) +
scale_fill_viridis_d() +
theme_minimal() +
guides(fill = FALSE) +
ylim(c(0, 1000)) +
labs(title = "Letters") 
Pregunta 4. Sesgos
Supón que los posibles sesgos pudieran estar representados por los predictores del dataset. ¿Podríamos determinar visual y numéricamente, con ayuda de PCA si los hay?
Para contrastar si nuestros datos son adecuados para hacer el analisis de componentes principales hay que pasar el test de esferecidad de Bartlett y el test de Kaiser-Meyer-Olkin.
En ambos tests rechazamos la hipotesis nula por lo que nuestros datos son adecuados para realizar el PCA Además el test KMO nos da un criterio de 0.68 que es un criterio adecuado para el PCA.
## Bartlett's Test of Sphericity
##
## Call: bart_spher(x = LetterRecognition[, -1])
##
## X2 = 166137.584
## df = 120
## p-value < 2.22e-16##
## Kaiser-Meyer-Olkin Statistics
##
## Call: KMOS(x = LetterRecognition[, -1])
##
## Measures of Sampling Adequacy (MSA):
## x.box y.box width high onpix x.bar y.bar x2bar
## 0.7268940 0.7790218 0.7430971 0.7566859 0.7086221 0.6109612 0.6396189 0.4546661
## y2bar xybar x2ybr xy2br x.ege xegvy y.ege yegvx
## 0.4716565 0.5183031 0.6176384 0.3643558 0.6728329 0.7434567 0.5081978 0.6030410
##
## KMO-Criterion: 0.6794291Vamos a ver la correlación que hay entre las variables predictoras para ver cómo se agrupan. Vemos que hay 2 grupos de variables muy correlados positivamente entre sí por lo que podriamos pensar que son candidatos a ser importantes en distintas componentes principales. Luego hay otro cluster más diverso en la correlación.
correlations <- cor(LetterRecognition[,c(-1)])
pheatmap(correlations,
clustering_method = "ward.D2",
color = colorRampPalette(c("#00AFBB", "white", "#FC4E07"))(255), # Paleta de colores
cellwidth = 12,
cellheight = 12,
breaks = seq(-1, 1, length.out = 256),
scale = "none",
annotation_legend = TRUE,
angle_col = 45, cutree_cols = 3)
Para aplicar el PCA de manera correcta, es esencial determinar si nuestros datos requieren un escalado apropiado.
datos_long_letter <- pivot_longer(LetterRecognition, cols = -lettr, names_to = "Variable", values_to = "Valor")
ggplot(datos_long_letter, aes(x = Variable, y = Valor, fill = Variable)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
scale_color_viridis(discrete = TRUE, alpha = 0.6) +
theme(
legend.position = "none",
plot.title = element_text(size = 14, hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 10),
axis.title.x = element_text(margin = margin(t = 10), size = 14, hjust = 0.5),
axis.title.y = element_text(margin = margin(r = 10), size = 14, hjust = 0.5)
) +
ggtitle("Variables") +
xlab("Variables") +
ylab("Value")
La distribución no es homogénea pero vemos que todos los datos están entre 0 y 15 por lo que no considero escalar los datos ya que todos los valores están restringidos entre esos valores y entiendo que tienen la misma escala.
var_exp <- pca_result.letter$sdev^2
prop_var_exp <- var_exp / sum(var_exp)
cum_var_exp <- cumsum(prop_var_exp)
df_var_exp <- data.frame(Comp = 1:length(prop_var_exp), VarExp = prop_var_exp)
df_cum_var_exp <- data.frame(Comp = 1:length(cum_var_exp), CumVarExp = cum_var_exp)
# Scree plot con los datos no escalados
ggplot(df_var_exp[1:10,], aes(x = Comp, y = VarExp)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_line(aes(group = 1), color = "blue") +
geom_point(color = "blue") +
theme_minimal() +
labs(x = "Principal components", y = "Variance", title = "Scree Plot") +
ylim(c(0,1)) +
geom_line(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), color="#8B1A1A") +
geom_point(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), color = "red") +
geom_bar(data = df_cum_var_exp[1:10,], aes(x = Comp, y = CumVarExp), stat = "identity", fill = "red", alpha= 0.25) +
annotate("text", x = 9, y = 0.75, label = "Cumulative Scree \n Plot", color = "#8B1A1A", size = 4) +
geom_text(data = df_cum_var_exp[1:10,], aes(x=Comp, y = CumVarExp +0.04, label = round(CumVarExp, 2)))
La primera componente solamente explica el 29% de la varianza. Con la segunda se explica hasta un 44% por lo que entendemos que las caracteristicas explicadas de nuestras variables están bien distribuidas y cada una explica algo importante en la clasificación de nuestras clases.
Para ver si hay sesgos claros en nuestros datos vamos a hacer una correlacion de las variables con los componentes principales.
sesgos.letter <- cor(pca_result.letter$x[,1:7], LetterRecognition[,-1])
pheatmap(sesgos.letter,
clustering_method = "ward.D2",
color = colorRampPalette(c("#00AFBB", "white", "#FC4E07"))(255), # Paleta de colores
cellwidth = 20,
cellheight = 20,
breaks = seq(-1, 1, length.out = 256),
scale = "none",
annotation_legend = TRUE,
cutree_cols = 3)
Vemos que las variables “x.ege”, “y.box”, “high”, “onpix”, “x.box” y “width” tienen una correlacion muy fuerte con la componente principal 1. Curiosamente, son las mismas variables que en el anterior mapa de correlación entre variables tenian una correlación muy fuerte entre ellas.
Vamos a ver la distribución de nuestra variables en los dos componentes principales 1 y 2.
pca.letter <- as.data.frame(pca_result.letter$x[,1:7], stringsAsFactors=F)
pca.letter <- cbind(LetterRecognition, pca.letter)
plots <- list()
for (i in 1:17) {
# Crea una gráfica para cada variable
p <- ggplot(pca.letter, aes_string(x = names(pca.letter)[18], y = names(pca.letter)[19], color = names(pca.letter)[i])) +
geom_point() +
labs(x = names(pca.letter)[18], y = names(pca.letter)[19], color = names(pca.letter)[i]) +
scale_color_gradient(low = "#D3D3D3", high = "#FC4E07") +
theme_minimal()
plots[[colnames(pca.letter)[i]]] <- p
}
plots$lettr <- ggplot(pca.letter, aes(x = PC1, y = PC2, color = lettr)) +
geom_point() +
labs(x = "PC 1", y = "PC 2", color = "lettr") +
theme_minimal()
grid.arrange(plots$x.box , plots$y.box, plots$width, plots$high, ncol=2)
Como vemos no se ven clusters separables en las dos primeras componentes principales pero sí vemos que cada variable se situa en unas zonas específicas, agrupandose entre ellas y con las variables más relacionadas entre ellas. Por ejemplo las variables “x.box”, “y.box”, “width” y “high” se encuentran en zonas parecidas ya que se correlacionan entre ellas.
Con respecto a la variable dependiente “lettr” vemos que podemos intuir que cada letra se separa de las otras en estas dos primeras componentes principales. En las siguientes ilustraciones muestro mejor la distribución de cada clase entre las dos primeras componentes principales.
unique_letters <- unique(LetterRecognition$lettr)
plots <- list()
for (letter in unique_letters) {
pca.letter$highlight <- ifelse(pca.letter$lettr == letter, "Highlighted", "Other")
plots[[letter]] <- ggplot(pca.letter, aes(x = PC1, y = PC2)) +
geom_point(aes(color = highlight), alpha = ifelse(pca.letter$highlight == "Highlighted", 1, 0.05)) +
scale_color_manual(values = c("Highlighted" = "#FC4E07", "Other" = "grey")) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 0, linetype = "dashed") +
labs(title = letter) +
theme_minimal() +
guides(color = FALSE, alpha = FALSE)
}
grid.arrange(plots$A, plots$B,plots$C,plots$D,plots$E,plots$F, plots$G, plots$H, plots$I,plots$J,ncol = 5)
grid.arrange(plots$K,plots$L, plots$M,plots$N,plots$O,plots$P,plots$Q,plots$R,plots$S,plots$T,ncol = 5)
Por último, otra forma de ver representada la distribución de las clases de letra pero esta vez entre las 3 primeras componentes principales. Se ve claramente que la clasificación de las clases de letras por medio de estas variables es muy prometedora.
plot <- plot_ly(data = pca.letter, x = ~PC1, y = ~PC2, z = ~PC3, color = ~lettr, type = 'scatter3d', mode = 'markers')
plot <- plot %>% layout(title = "3D PCA Letter Recongnition",
scene = list(xaxis = list(title = "PC1"),
yaxis = list(title = "PC2"),
zaxis = list(title = "PC3")))
plotLa gráfica anterior en 3 dimensiones puede ser explicada por la siguiente en la que mostramos el peso que tiene cada variable en cada dimensión. Por lo que “y.box” tiene el mayor peso a la hora de explicar la varianza de los datos en la componente principal 1. Las variables “x2ybr” y “y.bar” lo son importantes para la componente principal 2. Y las variables “x2bar” y “y2bar” para la componente principal 3.
p12 <- fviz_pca_var(pca_result.letter, col.var="cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, axes = c(1,2), label = "var", title="PC1-2") + theme(legend.position = "none")
p23 <- fviz_pca_var(pca_result.letter, col.var="cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, axes = c(2,3), label = "var", title="PC2-3") + theme(legend.position = "none")
grid.arrange(p12, p23, ncol= 2)
Ahora mostramos los valores numéricos de los vectores de carga numérica para las tres primeras componentes principales.
pca_rotation <- round(data.frame(pca_result.letter$rotation[,1:3]),2)
pca_rotation <- pca_rotation[order(-pca_rotation[,1]),]
formattable(pca_rotation,
list(
area(col= PC1:PC3) ~ formatter("span", style = x ~ ifelse(x > 0,
style(color = "green", font.weight = "bold"), style(color = "red", font.weight = "bold")))
))| PC1 | PC2 | PC3 | |
|---|---|---|---|
| y.box | 0.60 | -0.07 | -0.13 |
| high | 0.39 | -0.03 | 0.01 |
| onpix | 0.36 | 0.04 | 0.11 |
| width | 0.35 | -0.07 | 0.01 |
| x.box | 0.33 | -0.09 | 0.02 |
| x.ege | 0.25 | 0.01 | 0.39 |
| y.ege | 0.23 | 0.23 | -0.12 |
| x.bar | 0.05 | 0.26 | -0.06 |
| xybar | 0.04 | -0.23 | -0.45 |
| y2bar | 0.03 | 0.02 | -0.45 |
| x2bar | 0.00 | 0.17 | 0.57 |
| xegvy | 0.00 | -0.28 | 0.06 |
| xy2br | -0.01 | 0.13 | 0.03 |
| yegvx | -0.01 | 0.19 | 0.04 |
| y.bar | -0.03 | -0.54 | 0.04 |
| x2ybr | -0.05 | -0.60 | 0.23 |
Las variables con mayores cargas en la primera componente principañ (“x.box”, “y.box”, “width”, “high”, “onpix”) parecen estar relacionadas con las dimensiones físicas de las cajas que encierran las letras y la cantidad de píxeles encendidos (onpix), lo que sugiere que la PC1 podría estar capturando información relacionada con el tamaño y la ocupación espacial general de las letras.
Las variables con las mayores cargas absolutas en la segunda componente son “y.bar” y “x2ybr”, ambas con valores negativos, y “x.bar” con un valor positivo significativo. Esto indica que PC2 podría estar capturando variaciones en la distribución espacial de los píxeles en las letras, posiblemente reflejando características como la simetría o distribución vertical/horizontal de los píxeles.
Las variables “x2bar” y “xybar” tienen cargas significativas en la componente 3, aunque en direcciones opuestas. “x2bar” tiene una carga positiva, y xybar tiene una carga negativa. Estas cargas sugieren que la componente 3 podría estar asociada con la dispersión o la alineación de los píxeles en las imágenes de las letras, lo que podría interpretarse como características que capturan la textura o la estructura interna de las letras.
Ahora vamos a ver la explicación de cada variable a la varianza de cada componente principal.
p1 <- factoextra::fviz_cos2(pca_result.letter, choice = "var", axes = 1, top= 5, fill = "#00AFBB", ggtheme = theme_minimal()) + labs(y = "Cos2", title = "Dim 1")
p2 <- factoextra::fviz_cos2(pca_result.letter, choice = "var", axes = 2, top= 5, fill = "#00AFBB", ggtheme = theme_minimal())+ labs(y = "Cos2", title = "Dim 2")
p3 <- factoextra::fviz_cos2(pca_result.letter, choice = "var", axes = 3, top= 5, fill = "#00AFBB", ggtheme = theme_minimal())+ labs(y = "Cos2", title = "Dim 3")
grid.arrange(p1, p2, p3, ncol=3)
Las cinco variables que mejor se proyectan en la primera componente principal son mostradas en el primer gráfico. La variable “y.box” tiene el valor de Cos2 más alto, lo que indica que es la que mejor se representa en esta dimensión. Esto significa que “y.box” contribuye significativamente a la varianza capturada por la primera componente principal.
El segundo gráfico muestra las cinco variables con el mayor Cos2 para la segunda componente principal. La variable “x2ybar” tiene un valor de Cos2 mucho mayor en comparación con las otras cuatro, lo que sugiere que tiene una influencia preponderante en esta segunda dimensión.
La tercera gráfica se presentan las cinco variables más importantes en términos de representación en la tercera componente principal. La variable “x2ybar” también aparece aquí con el valor más alto, seguida de ‘y2bar’, ‘xybar’, ‘x.ege’, y ‘x2ybr’, lo que sugiere que estas variables también son importantes en la definición de esta tercera componente principal.
Pregunta 5. Elección algoritmo
¿Qué algoritmo de entre los que conoces aplicarías según el problema para el que lo has identificado?
KNN es una elección razonable para el conjunto de datos letterRecognition debido a su capacidad para manejar clasificaciones multiclase y su utilidad en situaciones donde la clasificación depende de la similitud de las instancias en un espacio de características. Su robustez frente a las variaciones complejas en los datos lo hace adecuado para reconocer patrones visuales, lo cual es esencial en la clasificación de letras basada en atributos visuales.
Además, KNN no hace suposiciones sobre la distribución subyacente de los datos, lo cual es ventajoso en el procesamiento de imágenes, donde las distribuciones pueden ser complicadas y no lineales. La capacidad de KNN para adaptarse a los datos localmente, ajustando el número de vecinos, permite una gran flexibilidad para encontrar estructuras en el espacio de características que son significativas para la clasificación.
Sin embargo, es importante destacar que KNN tiene sus desventajas. La principal es su naturaleza computacionalmente intensiva, especialmente en conjuntos de datos grandes con muchas características, debido a la necesidad de calcular la distancia entre pares para cada instancia de prueba contra todas las instancias de entrenamiento. Además, KNN puede ser sensible a características irrelevantes o redundantes, por lo que un buen preprocesamiento y posiblemente la selección de características son pasos cruciales para su rendimiento óptimo.
Pregunta 6. Transformación
¿Sería necesario aplicar alguna transformación a los datos? ¿Por qué? Si la respuesta es afirmativa, explica qué tendría de beneficioso hacerlo.
En este caso, los datos fueron escalados previamente para situarlos en un rango entre 0 y 15, por lo que no consideraría escalar los datos.
Por otro lado, podría reducir la dimensionalidad a 7 u 8 dimensiones ya que explican más del 80% de la varianza de los datos en el caso de que la capacidad computacional se viese comprometida. En este caso, no voy a reducir la dimensionalidad y voy a seguir con los datos sin transformar.
Pregunta 7. Algoritmo
Aplica el algoritmo designado en la pregunta 5 y optimiza el estadístico o estadísticos que consideres pertinente para medir su rendimiento mediante la selección de un modelo óptimo. Comenta los resultados
En primer lugar, vamos a generar los datos de entrenamiento y test de los datos de letterRecognition para entrenar el algoritmo de KNN. Vemos que el porcentaje de cada letra en entrenamiento y test se mantiene.
set.seed(1234)
letterRecognition.TrainIdx.80<- createDataPartition(LetterRecognition$lettr,
p=0.8, #Genera un 80% para train, 20% para test
list = FALSE, #Dame los resultados en una matriz
times = 1) #Genera solamente una partición 80/20
trainSet.letter <- LetterRecognition[letterRecognition.TrainIdx.80, ]
testSet.letter <- LetterRecognition[-letterRecognition.TrainIdx.80, ]
table(trainSet.letter$lettr)##
## A B C D E F G H I J K L M N O P Q R S T
## 632 613 589 644 615 620 619 588 604 598 592 609 634 627 603 643 627 607 599 637
## U V W X Y Z
## 651 612 602 630 629 588##
## A B C D E F G H I J K L M N O P Q R S T
## 157 153 147 161 153 155 154 146 151 149 147 152 158 156 150 160 156 151 149 159
## U V W X Y Z
## 162 152 150 157 157 146En el apartado anterior de los datos BostonHousing explico la información del paquete kknn de caret y los hiperparámetros a utilizar (kmax, distance y kernel).
Para este modelo he elegido entrenar los datos con validación cruzada sin repetición, ya que la dimensionalidad de estos datos es bastante grande y a nivel computacional es más costoso. Permito el análisis en paralelo para aumentar el rendimiento
set.seed(1234)
kknnControl.letter<- trainControl(method = "cv",
number = 10 ,
returnResamp = "final",
seeds = NULL,
allowParallel = T)En un primer intento de buscar el mejor modelo realizo una parrilla de hiperparámetros con los siguientes valores dando 9 posibilidades combinatorias:
kmax: 5, 7, 9
distance: 1, 2, 3
kernel: optimal
set.seed(1234)
y <- trainSet.letter$lettr
x <- trainSet.letter[,-1]
mygrid.letter <- expand.grid(kmax = seq(5,9,2),
distance = c(1,2,3),
kernel = c(#"rectangular","triangular",
"optimal"
# ,"epanechnikov", "biweight", "triweight"
# ,"cos", "inv", "gaussian","rank"
))
myfit.letter.cv10 <- train(y=y, x=x,
method = "kknn",
metric = "Kappa",
trControl = kknnControl.letter,
tuneGrid=mygrid.letter
)
(myfit.letter.cv10.best <- subset(myfit.letter.cv10$results, kmax == myfit.letter.cv10$bestTune$kmax & distance == myfit.letter.cv10$bestTune$distance & kernel == myfit.letter.cv10$bestTune$kernel))
En este primer modelo vemos que los mejores hiperparámetros son con una kmax de 9, distancia Manhattan de 1 y el kernel optimal. La Accuracy (0.95) y Kappa (0.95) dan valores muy buenos para este modelo.
Vamos a probar con kmax más altos ya que parece que la tendencia del kmax es a subir. Tambien probar otros tipos de kernel para ver si mejoramos el modelo. Los hiperparámetros son:
kmax: 9, 10, 11
distance: 1
kernel: rectangular, triangular, optimal
mygrid.letter <- expand.grid(kmax = seq(9,11,1),
distance = 1,
kernel = c("rectangular","triangular",
"optimal"
# ,"epanechnikov", "biweight", "triweight"
# ,"cos", "inv", "gaussian","rank"
))
myfit.letter.cv10.2 <- train(y=y, x=x,
method = "kknn",
metric = "Kappa",
trControl = kknnControl.letter,
tuneGrid=mygrid.letter
)
(myfit.letter.cv10.2.best <- subset(myfit.letter.cv10.2$results, kmax == myfit.letter.cv10.2$bestTune$kmax & distance == myfit.letter.cv10.2$bestTune$distance & kernel == myfit.letter.cv10.2$bestTune$kernel))
En este segundo modelo vemos que los mejores hiperparámetros sigue siendo una kmax de 9 y el kernel triangular. La Accuracy (0.96) y Kappa (0.96) dan valores muy buenos para este modelo, mejorando los valores anteriores.
Vamos a probar con kmax estable en 9 pero con otros tipos de kernel para ver si mejoramos el modelo. Los hiperparámetros son:
kmax: 9,
distance: 1
kernel: triangular, epanechnikov, biweight, triweight, cos, inv, gaussian, rank
mygrid.letter <- expand.grid(kmax = 9,
distance = 1,
kernel = c(#"rectangular",
"triangular"
#,"optimal"
,"epanechnikov", "biweight", "triweight"
,"cos", "inv", "gaussian","rank"
))
myfit.letter.cv10.3 <- train(y=y, x=x,
method = "kknn",
metric = "Kappa",
trControl = kknnControl.letter,
tuneGrid=mygrid.letter
)
(myfit.letter.cv10.3.best <- subset(myfit.letter.cv10.3$results, kmax == myfit.letter.cv10.3$bestTune$kmax & distance == myfit.letter.cv10.3$bestTune$distance & kernel == myfit.letter.cv10.3$bestTune$kernel))
Seguimos viendo que otros kernel no mejoramos el modelo así que el mejor modelo para este dataset es con una kmax de 9, distancia Manhattan de 1 y kernel triangular.
Como vemos en las siguiente gráficas, casi todas las iteraciones de nuestro mejor modelo de entrenamiento nos da como resultado un accuracy de entorno 0.96 y un Kappa de 0.96. Estos valores son muy altos y muy buenos, por lo que deberiamos mirar si nuestro modelo ha caido en overfitting.
resample <- myfit.letter.cv10.3$resample
p1 <- ggplot(data = resample, aes(x = Accuracy)) +
geom_density(size = 2,color = "#8B1A1A") +
theme_minimal() +
labs(title = "Accuracy",
x = "Accuracy")+
geom_vline(xintercept = median(resample$Accuracy),
linetype = "dashed") +
annotate("text", x = 0.95, y = 4,
label = paste("Accuracy", "\n", round(median(resample$Accuracy), 3)), color = "#8B1A1A", size = 4) +
xlim(c(0.93, 0.97))
p2 <- ggplot(data = resample, aes(x = Kappa)) +
geom_density(size = 2, color = "#f1e6b2") +
theme_minimal() +
labs(title = "Kappa",
x = "Kappa")+
geom_vline(xintercept = median(resample$Kappa),
linetype = "dashed") +
annotate("text", x = 0.95, y = 6,
label = paste("Kappa", "\n", round(median(resample$Kappa), 3)), color = "#c0b88e", size = 4) +
xlim(c(0.93, 0.97))
grid.arrange(p1, p2, ncol = 2)
En la siguiente tabla vemos las variables más importantes para nuestro modelo para cada una de las letras de nuestra variable dependiente lettr. La variable “y2bar” tiene un fuerte importancia en cada una de las letras. En cambio, “y.box” no es importante en nuestro modelo para predecir letras a partir de estas variables.
import <- varImp(myfit.letter.cv10.3)
import$importance <- round(import$importance[order(-import$importance[,1]),], 2)
formattable(import$importance, list(
area(col= A:Z) ~ color_bar("lightblue")
))| A | B | C | D | E | F | G | H | I | J | K | L | M | N | O | P | Q | R | S | T | U | V | W | X | Y | Z | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| y2bar | 97.52 | 99.70 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 98.84 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 99.33 | 96.08 | 96.08 | 96.08 | 96.08 | 96.08 | 97.10 | 96.08 | 99.70 |
| x2ybr | 91.79 | 79.80 | 84.54 | 92.69 | 84.14 | 80.88 | 81.44 | 79.80 | 81.67 | 79.80 | 84.79 | 79.80 | 82.39 | 82.40 | 85.35 | 79.80 | 79.80 | 98.03 | 97.78 | 99.92 | 96.75 | 80.82 | 100.00 | 81.49 | 79.80 | 91.79 |
| x.bar | 90.05 | 36.65 | 69.71 | 80.71 | 62.23 | 45.86 | 46.29 | 35.22 | 76.66 | 78.37 | 35.22 | 54.79 | 46.37 | 68.17 | 35.22 | 48.11 | 35.22 | 69.25 | 69.26 | 65.04 | 63.72 | 48.25 | 55.17 | 43.14 | 82.49 | 90.05 |
| x2bar | 89.67 | 84.70 | 76.48 | 76.48 | 91.44 | 76.48 | 76.48 | 76.48 | 76.48 | 76.48 | 76.48 | 76.48 | 96.14 | 78.06 | 91.34 | 77.63 | 76.48 | 76.48 | 92.65 | 76.48 | 76.48 | 76.48 | 76.48 | 76.48 | 76.48 | 89.67 |
| xegvy | 89.62 | 75.50 | 87.69 | 93.32 | 81.41 | 80.29 | 85.64 | 75.50 | 75.50 | 75.50 | 75.50 | 86.69 | 84.36 | 97.22 | 77.25 | 75.50 | 75.98 | 92.12 | 88.18 | 98.22 | 97.53 | 79.52 | 99.31 | 79.93 | 75.50 | 89.62 |
| y.ege | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 89.61 | 93.39 | 89.61 | 53.53 |
| y.bar | 89.40 | 84.08 | 87.04 | 99.78 | 84.08 | 85.86 | 84.08 | 84.08 | 84.27 | 84.08 | 84.08 | 91.27 | 84.08 | 98.63 | 84.08 | 90.97 | 84.08 | 99.49 | 84.08 | 98.23 | 96.11 | 85.51 | 96.56 | 84.08 | 84.08 | 89.40 |
| xy2br | 86.91 | 89.95 | 79.25 | 83.63 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 95.59 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 79.25 | 99.12 | 89.95 |
| yegvx | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 65.63 | 42.28 |
| x.ege | 42.96 | 32.26 | 18.04 | 46.46 | 8.17 | 53.03 | 80.15 | 69.74 | 57.49 | 66.70 | 98.82 | 86.01 | 29.41 | 35.11 | 27.54 | 36.65 | 34.46 | 46.61 | 50.79 | 20.58 | 97.76 | 21.47 | 36.98 | 64.77 | 46.43 | 42.96 |
| onpix | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 49.78 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 42.68 | 45.54 | 42.68 | 42.68 | 42.68 | 49.65 | 31.12 |
| xybar | 25.45 | 9.66 | 16.81 | 83.43 | 13.81 | 10.74 | 27.22 | 84.02 | 12.84 | 44.15 | 8.82 | 11.49 | 6.92 | 68.76 | 34.53 | 3.64 | 9.16 | 22.38 | 10.05 | 3.64 | 4.47 | 3.64 | 3.64 | 38.23 | 22.82 | 25.45 |
| x.box | 24.95 | 26.80 | 23.29 | 23.29 | 25.45 | 34.49 | 40.85 | 23.29 | 38.43 | 23.29 | 47.71 | 36.93 | 28.18 | 28.55 | 25.42 | 28.26 | 27.76 | 23.29 | 39.74 | 28.30 | 54.59 | 25.03 | 25.05 | 23.29 | 23.29 | 26.80 |
| width | 13.93 | 2.92 | 14.56 | 7.98 | 7.62 | 15.75 | 62.68 | 39.05 | 21.69 | 28.30 | 33.50 | 16.68 | 9.44 | 3.37 | 7.37 | 3.70 | 3.77 | 11.59 | 5.37 | 4.49 | 41.93 | 12.31 | 4.82 | 10.90 | 9.48 | 13.93 |
| high | 3.96 | 3.61 | 1.09 | 2.08 | 1.16 | 1.09 | 2.62 | 6.24 | 4.45 | 3.12 | 1.89 | 4.88 | 4.02 | 8.00 | 27.84 | 1.86 | 7.01 | 1.09 | 1.09 | 3.71 | 5.88 | 1.96 | 11.54 | 7.63 | 2.38 | 3.96 |
| y.box | 2.14 | 3.50 | 0.00 | 0.44 | 1.16 | 1.56 | 0.11 | 4.91 | 4.09 | 1.29 | 1.64 | 3.16 | 2.05 | 2.18 | 14.09 | 1.54 | 7.25 | 0.88 | 0.08 | 3.99 | 4.19 | 0.77 | 5.53 | 5.53 | 1.80 | 3.50 |
Vamos ahora a predecir las clases de los datos que tenemos de test para este dataset que hemos generado anteriormente. Despues de eso vamos mostrar la información de la matriz de confusión de estos resultados.
Una matriz de confusión es una herramienta poderosa utilizada en el campo del aprendizaje automático para evaluar el rendimiento de los modelos de clasificación. Es una tabla específica que permite la visualización del desempeño de un algoritmo de clasificación, mostrando de manera explícita cuándo las clases son confundidas por el modelo durante la predicción.
La matriz de confusión compara las etiquetas predichas por el modelo contra las etiquetas reales (verdaderas) que se encuentran en el conjunto de datos de prueba. La matriz de confusión permite calcular métricas importantes como la precisión, la sensibilidad, la especificidad, y el valor predictivo. Tambien ayuda a entender los tipos de errores que está cometiendo el modelo (por ejemplo, si está prediciendo falsos positivos más a menudo que falsos negativos).
preds.letter <- predict(myfit.letter.cv10.3$finalModel, newdata = testSet.letter)
conf.letter <- confusionMatrix(preds.letter, testSet.letter$lettr)
conf.letter## Confusion Matrix and Statistics
##
## Reference
## Prediction A B C D E F G H I J K L M N O P Q
## A 155 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## B 0 145 0 2 0 1 0 2 0 0 0 0 3 3 1 1 0
## C 0 0 144 0 1 0 0 0 0 0 0 0 0 0 1 0 0
## D 0 0 0 156 0 0 3 3 0 1 1 0 1 0 4 0 0
## E 0 1 0 0 144 1 3 0 1 0 0 0 0 0 0 0 0
## F 0 0 0 0 0 147 0 0 0 0 0 0 0 0 0 9 0
## G 0 0 0 0 1 0 147 2 0 0 0 0 0 0 1 0 0
## H 0 2 0 1 1 0 0 128 0 0 2 1 0 0 0 0 0
## I 0 0 0 0 0 1 0 0 144 7 0 0 0 0 0 0 0
## J 0 0 0 0 0 0 0 0 6 141 0 0 0 0 0 0 0
## K 0 0 0 0 0 0 0 3 0 0 135 0 0 0 0 0 0
## L 0 0 1 0 1 0 0 0 0 0 0 149 0 0 0 1 0
## M 2 0 0 0 0 0 0 1 0 0 0 0 148 1 0 0 0
## N 0 0 0 2 0 1 0 0 0 0 0 0 1 147 0 0 0
## O 0 0 1 0 0 0 0 2 0 0 0 0 1 1 142 0 8
## P 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 148 0
## Q 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 146
## R 0 2 0 0 0 0 0 4 0 0 2 0 0 1 0 0 1
## S 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## T 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## U 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0
## V 0 2 0 0 0 0 0 0 0 0 0 0 2 3 0 0 0
## W 0 0 1 0 0 0 0 0 0 0 0 0 2 0 1 0 0
## X 0 0 0 0 2 0 1 0 0 0 5 1 0 0 0 1 0
## Y 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1
## Z 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0
## Reference
## Prediction R S T U V W X Y Z
## A 0 0 0 0 0 0 0 1 0
## B 4 0 0 0 1 0 0 0 0
## C 0 0 0 0 0 0 0 0 0
## D 0 0 0 0 0 0 1 0 0
## E 0 0 0 0 0 0 0 0 0
## F 0 0 0 0 0 0 0 0 0
## G 1 0 0 0 0 0 0 0 0
## H 1 0 0 0 0 0 0 0 0
## I 0 0 0 0 0 0 0 0 0
## J 0 0 0 0 0 0 0 0 0
## K 4 0 0 0 0 0 1 0 0
## L 0 0 0 0 0 0 0 0 0
## M 0 0 0 0 1 0 0 0 0
## N 2 0 0 0 0 1 0 0 0
## O 0 0 0 0 1 0 0 0 0
## P 0 0 0 0 0 0 0 1 0
## Q 1 0 0 0 0 0 0 1 2
## R 138 0 0 0 0 0 3 0 0
## S 0 149 0 0 0 0 0 0 0
## T 0 0 156 0 0 0 0 1 0
## U 0 0 0 162 0 1 0 0 0
## V 0 0 0 0 148 0 0 1 0
## W 0 0 0 0 0 148 0 0 0
## X 0 0 0 0 0 0 152 1 0
## Y 0 0 3 0 1 0 0 151 0
## Z 0 0 0 0 0 0 0 0 144
##
## Overall Statistics
##
## Accuracy : 0.9564
## 95% CI : (0.9496, 0.9625)
## No Information Rate : 0.0406
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9546
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: A Class: B Class: C Class: D Class: E Class: F
## Sensitivity 0.98726 0.94771 0.97959 0.96894 0.94118 0.94839
## Specificity 0.99974 0.99531 0.99948 0.99634 0.99844 0.99765
## Pos Pred Value 0.99359 0.88957 0.98630 0.91765 0.96000 0.94231
## Neg Pred Value 0.99948 0.99791 0.99922 0.99869 0.99766 0.99791
## Prevalence 0.03937 0.03837 0.03686 0.04037 0.03837 0.03887
## Detection Rate 0.03887 0.03636 0.03611 0.03912 0.03611 0.03686
## Detection Prevalence 0.03912 0.04087 0.03661 0.04263 0.03761 0.03912
## Balanced Accuracy 0.99350 0.97151 0.98954 0.98264 0.96981 0.97302
## Class: G Class: H Class: I Class: J Class: K Class: L
## Sensitivity 0.95455 0.87671 0.95364 0.94631 0.91837 0.98026
## Specificity 0.99870 0.99792 0.99792 0.99844 0.99792 0.99922
## Pos Pred Value 0.96711 0.94118 0.94737 0.95918 0.94406 0.98026
## Neg Pred Value 0.99818 0.99533 0.99818 0.99792 0.99688 0.99922
## Prevalence 0.03862 0.03661 0.03786 0.03736 0.03686 0.03811
## Detection Rate 0.03686 0.03210 0.03611 0.03536 0.03385 0.03736
## Detection Prevalence 0.03811 0.03410 0.03811 0.03686 0.03586 0.03811
## Balanced Accuracy 0.97662 0.93732 0.97578 0.97237 0.95814 0.98974
## Class: M Class: N Class: O Class: P Class: Q Class: R
## Sensitivity 0.93671 0.94231 0.94667 0.92500 0.93590 0.91391
## Specificity 0.99869 0.99817 0.99635 0.99896 0.99870 0.99661
## Pos Pred Value 0.96732 0.95455 0.91026 0.97368 0.96689 0.91391
## Neg Pred Value 0.99739 0.99765 0.99791 0.99687 0.99739 0.99661
## Prevalence 0.03962 0.03912 0.03761 0.04012 0.03912 0.03786
## Detection Rate 0.03711 0.03686 0.03561 0.03711 0.03661 0.03460
## Detection Prevalence 0.03837 0.03862 0.03912 0.03811 0.03786 0.03786
## Balanced Accuracy 0.96770 0.97024 0.97151 0.96198 0.96730 0.95526
## Class: S Class: T Class: U Class: V Class: W Class: X
## Sensitivity 1.00000 0.98113 1.00000 0.97368 0.98667 0.96815
## Specificity 0.99974 0.99948 0.99922 0.99791 0.99896 0.99713
## Pos Pred Value 0.99333 0.98734 0.98182 0.94872 0.97368 0.93252
## Neg Pred Value 1.00000 0.99922 1.00000 0.99896 0.99948 0.99869
## Prevalence 0.03736 0.03987 0.04062 0.03811 0.03761 0.03937
## Detection Rate 0.03736 0.03912 0.04062 0.03711 0.03711 0.03811
## Detection Prevalence 0.03761 0.03962 0.04137 0.03912 0.03811 0.04087
## Balanced Accuracy 0.99987 0.99030 0.99961 0.98580 0.99281 0.98264
## Class: Y Class: Z
## Sensitivity 0.96178 0.98630
## Specificity 0.99843 0.99922
## Pos Pred Value 0.96178 0.97959
## Neg Pred Value 0.99843 0.99948
## Prevalence 0.03937 0.03661
## Detection Rate 0.03786 0.03611
## Detection Prevalence 0.03937 0.03686
## Balanced Accuracy 0.98011 0.99276La accuracy global del modelo es del 95.64%, indicando un alto nivel de precisión en la clasificación de las instancias en sus respectivas categorías. El valor de Kappa es de 0.9546, lo que indica una muy buena concordancia entre las predicciones del modelo y los valores reales, más allá de lo que se esperaría por casualidad.
Tambien cada clase muestra valores altos de sensibilidad (la capacidad del modelo para identificar correctamente los verdaderos positivos) y especificidad (la capacidad del modelo para identificar correctamente los verdaderos negativos), lo que demuestra que el modelo es capaz de clasificar correctamente la mayoría de las instancias en todas las categorías.
Aunque el modelo es super robusto y comete muy pocos errores, vamos a ver, mediante el siguiente heatmap, cuáles son los fallos que comete el modelo.
conf_matrix_long <- as.data.frame(conf.letter$table)
ggplot(data = conf_matrix_long, aes(x = Reference, y = Prediction, fill = Freq)) +
geom_tile() +
scale_fill_gradient2(low = "white", high = "#FC4E07", mid ="#00AFBB", midpoint = 5, limit = c(0, 10)) +
theme_minimal() +
labs(x = "Predicted", y = "Reference", fill = "Count") +
geom_text(aes(x=Reference, y = Prediction, label=ifelse(Freq > 0, paste(Reference, Prediction), "")), size = 2)+ coord_flip()
Lo que vemos en el heatmap son los fallos cometimos de nuestro modelo. Por ejemplo, ha cometido casi 10 fallos al predecir una letra P en lugar de una F. Otro fallo que ha cometido alguna que otra vez es predecir una Q en lugar de una O. Son fallos por una similitud alta a la hora de manuscribir las letras.
Aun con estos pocos fallos, nuestro modelo predice perfectamente con este modelo de entrenamiento las distintas clases. Vamos a ver a continuación si este modelo ha incurrido en overfitting.
# Evaluar el rendimiento en el conjunto de entrenamiento
trainPredictions <- predict(myfit.letter.cv10.3$finalModel, trainSet.letter)
conf_train <- confusionMatrix(trainPredictions, trainSet.letter$lettr)
# Evaluar el rendimiento en el conjunto de prueba
testPredictions <- predict(myfit.letter.cv10.3$finalModel, testSet.letter)
conf_test <- confusionMatrix(testPredictions, testSet.letter$lettr)
conf_train$overall## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.99862603 0.99857102 0.99792053 0.99913875 0.04065701
## AccuracyPValue McnemarPValue
## 0.00000000 NaN## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.95636911 0.95462200 0.94956105 0.96249780 0.04062187
## AccuracyPValue McnemarPValue
## 0.00000000 NaNAl calcular el accuracy de los valores predichos con los datos de entrenamiento y el accuracy de los datos de test, observamos que en los datos de entrenamiento la precisión es mayor que en los datos de test pero la diferencia no es muy grande, siendo el accuracy también muy alto. Encontramos un modelo muy bueno para predecir el tipo de letra según las variables recogidas en este dataset.
Diabetes BRFSS
Los datos de este dataset se localizan en el archivo csv “diabetesBRFSS2015.csv”.
Pregunta 1. Tamaño y tipo
¿Qué tamaño tiene? ¿De qué tipo son las variables?
## 'data.frame': 70692 obs. of 22 variables:
## $ Diabetes_binary : num 0 0 0 0 0 0 0 0 0 0 ...
## $ HighBP : num 1 1 0 1 0 0 0 0 0 0 ...
## $ HighChol : num 0 1 0 1 0 0 1 0 0 0 ...
## $ CholCheck : num 1 1 1 1 1 1 1 1 1 1 ...
## $ BMI : num 26 26 26 28 29 18 26 31 32 27 ...
## $ Smoker : num 0 1 0 1 1 0 1 1 0 1 ...
## $ Stroke : num 0 1 0 0 0 0 0 0 0 0 ...
## $ HeartDiseaseorAttack: num 0 0 0 0 0 0 0 0 0 0 ...
## $ PhysActivity : num 1 0 1 1 1 1 1 0 1 0 ...
## $ Fruits : num 0 1 1 1 1 1 1 1 1 1 ...
## $ Veggies : num 1 0 1 1 1 1 1 1 1 1 ...
## $ HvyAlcoholConsump : num 0 0 0 0 0 0 1 0 0 0 ...
## $ AnyHealthcare : num 1 1 1 1 1 0 1 1 1 1 ...
## $ NoDocbcCost : num 0 0 0 0 0 0 0 0 0 0 ...
## $ GenHlth : num 3 3 1 3 2 2 1 4 3 3 ...
## $ MentHlth : num 5 0 0 0 0 7 0 0 0 0 ...
## $ PhysHlth : num 30 0 10 3 0 0 0 0 0 6 ...
## $ DiffWalk : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Sex : num 1 1 1 1 0 0 1 1 0 1 ...
## $ Age : num 4 12 13 11 8 1 13 6 3 6 ...
## $ Education : num 6 6 6 6 5 4 5 4 6 4 ...
## $ Income : num 8 8 8 8 8 7 6 3 8 4 ...Es un dataset con 70,692 muestras y 22 variables casi todas numéricas aunque es un problema de la carga de los datos. Paso las variables categoricas a factorizarlas.
diabetesBRFSS2015[] <- lapply(diabetesBRFSS2015, as.factor)
diabetesBRFSS2015$BMI <- as.numeric(diabetesBRFSS2015$BMI)
diabetesBRFSS2015$MentHlth <- as.numeric(diabetesBRFSS2015$MentHlth)
diabetesBRFSS2015$PhysHlth <- as.numeric(diabetesBRFSS2015$PhysHlth)
diabetesBRFSS2015$Diabetes_binary <- factor(diabetesBRFSS2015$Diabetes_binary, levels= c(0,1), labels = c("Healthy", "Diabetic"))
Hmisc::describe(diabetesBRFSS2015)## diabetesBRFSS2015
##
## 22 Variables 70692 Observations
## --------------------------------------------------------------------------------
## Diabetes_binary
## n missing distinct
## 70692 0 2
##
## Value Healthy Diabetic
## Frequency 35346 35346
## Proportion 0.5 0.5
## --------------------------------------------------------------------------------
## HighBP
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 30860 39832
## Proportion 0.437 0.563
## --------------------------------------------------------------------------------
## HighChol
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 33529 37163
## Proportion 0.474 0.526
## --------------------------------------------------------------------------------
## CholCheck
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 1749 68943
## Proportion 0.025 0.975
## --------------------------------------------------------------------------------
## BMI
## n missing distinct Info Mean Gmd .05 .10
## 70692 0 80 0.997 18.86 7.422 10 11
## .25 .50 .75 .90 .95
## 14 18 22 28 32
##
## lowest : 1 2 3 4 5, highest: 76 77 78 79 80
## --------------------------------------------------------------------------------
## Smoker
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 37094 33598
## Proportion 0.525 0.475
## --------------------------------------------------------------------------------
## Stroke
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 66297 4395
## Proportion 0.938 0.062
## --------------------------------------------------------------------------------
## HeartDiseaseorAttack
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 60243 10449
## Proportion 0.852 0.148
## --------------------------------------------------------------------------------
## PhysActivity
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 20993 49699
## Proportion 0.297 0.703
## --------------------------------------------------------------------------------
## Fruits
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 27443 43249
## Proportion 0.388 0.612
## --------------------------------------------------------------------------------
## Veggies
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 14932 55760
## Proportion 0.211 0.789
## --------------------------------------------------------------------------------
## HvyAlcoholConsump
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 67672 3020
## Proportion 0.957 0.043
## --------------------------------------------------------------------------------
## AnyHealthcare
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 3184 67508
## Proportion 0.045 0.955
## --------------------------------------------------------------------------------
## NoDocbcCost
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 64053 6639
## Proportion 0.906 0.094
## --------------------------------------------------------------------------------
## GenHlth
## n missing distinct
## 70692 0 5
##
## Value 1 2 3 4 5
## Frequency 8282 19872 23427 13303 5808
## Proportion 0.117 0.281 0.331 0.188 0.082
## --------------------------------------------------------------------------------
## MentHlth
## n missing distinct Info Mean Gmd .05 .10
## 70692 0 31 0.685 4.752 6.285 1 1
## .25 .50 .75 .90 .95
## 1 1 3 16 31
##
## lowest : 1 2 3 4 5, highest: 27 28 29 30 31
## --------------------------------------------------------------------------------
## PhysHlth
## n missing distinct Info Mean Gmd .05 .10
## 70692 0 31 0.818 6.81 8.949 1 1
## .25 .50 .75 .90 .95
## 1 1 7 31 31
##
## lowest : 1 2 3 4 5, highest: 27 28 29 30 31
## --------------------------------------------------------------------------------
## DiffWalk
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 52826 17866
## Proportion 0.747 0.253
## --------------------------------------------------------------------------------
## Sex
## n missing distinct
## 70692 0 2
##
## Value 0 1
## Frequency 38386 32306
## Proportion 0.543 0.457
## --------------------------------------------------------------------------------
## Age
## n missing distinct
## 70692 0 13
##
## Value 1 2 3 4 5 6 7 8 9 10 11
## Frequency 979 1396 2049 2793 3520 4648 6872 8603 10112 10856 8044
## Proportion 0.014 0.020 0.029 0.040 0.050 0.066 0.097 0.122 0.143 0.154 0.114
##
## Value 12 13
## Frequency 5394 5426
## Proportion 0.076 0.077
## --------------------------------------------------------------------------------
## Education
## n missing distinct
## 70692 0 6
##
## Value 1 2 3 4 5 6
## Frequency 75 1647 3447 19473 20030 26020
## Proportion 0.001 0.023 0.049 0.275 0.283 0.368
## --------------------------------------------------------------------------------
## Income
## n missing distinct
## 70692 0 8
##
## Value 1 2 3 4 5 6 7 8
## Frequency 3611 4498 5557 6658 8010 10287 11425 20646
## Proportion 0.051 0.064 0.079 0.094 0.113 0.146 0.162 0.292
## --------------------------------------------------------------------------------No hay datos nulos en nuestro dataset. Hay algunas variables categóricas bastante balanceadas pero otras no. Parece que hay algunas variables con varias clases e incluso otras que pueden ser numéricas.
Pregunta 2. Representación ejemplos
Explica qué representan los ejemplos
Los datos del dataset de diabetesBRFSS2015 representa información de encuestas sobre estilos de vida de personas en general, incluyendo su diagnóstico de diabetes. Cada instancia en este conjunto de datos representa a una persona que participó en este estudio. El objetivo principal de este conjunto de datos es entender mejor la relación entre el estilo de vida y la diabetes en EE.UU.
El conjunto de datos incluye 22 variables que abarcan desde información demográfica hasta resultados de exámenes de laboratorio y respuestas a preguntas de encuestas. Aquí hay un resumen de algunas de las variables clave:
- Diabetes_binary: Variable objetivo binaria que indica si el paciente tiene diabetes (incluyendo prediabetes) o no.
- HighBP: Presión arterial alta, indicada binariamente 0 = no 1 = si.
- HighChol: Colesterol alto, también indicado de manera binaria 0 = no 1 = si.
- CholCheck: Indica si el paciente ha tenido un chequeo de colesterol en los últimos 5 años 0 = no 1 = si.
- BMI: Índice de masa corporal, medida cuantitativa.
- Smoker: Indica si el paciente ha fumado al menos 100 cigarrillos en su vida 0 = no 1 = si.
- Stroke: Indica si al paciente le han diagnosticado un derrame cerebral 0 = no 1 = si.
- HeartDiseaseorAttack: Presencia de enfermedad coronaria o ataque cardíaco, indicado binariamente 0 = no 1 = si.
- PhysActivity: Actividad física en los últimos 30 días, excluyendo el trabajo 0 = no 1 = si, indicado binariamente.
- Fruits: Consumo diario de frutas, indicado binariamente 0 = no 1 = si.
- Veggies: Consumo diario verduras, indicado binariamente 0 = no 1 = si.
- HvyAlcoholConsump: Consumo excesivo de alcohol 0 = no 1 = si.
- AnyHealthcare: Cobertura de salud de cualquier tipo 0 = no 1 = si.
- NoDocbcCost: Si hubo un momento en los últimos 12 meses en que el paciente necesitaba ver a un doctor pero no pudo debido al costo 0 = no 1 = si.
- GenHlth: Estado general de salud, medido en una escala de 1 a 5, 1 = excelente 2 = muy buena 3 = buena 4 = poca 5 = muy poca.
- MentHlth: Salud mental indicada por el número de días en los últimos 30 días que no fueron buenos.
- PhysHlth: Salud física indicada por el número de días en los últimos 30 días que no fueron buenos.
- DiffWalk: Dificultad para caminar o subir escaleras. 0 = no 1 = si.
- Sex: Género del paciente. 0 = mujer 1 = hombre
- Age: Edad, categorizada en 13 niveles de edad. 1 = 18-24, 9 = 60-64, 13 = 80 o mayores.
- Education: Nivel de educación, en una escala de 1 a 6, 1 = Nunca fue a la escuela, 2 = Grados 1 a 8 (Elemental), 3 = Grados 9 a 11 (Algo de instituto) 4 = Grados 12 or graduación (Graduado de instituto) 5 = Universidad de 1 a 3 años (Asistió a la universidad) 6 = Universidad 4 años o más (Graduado universitario).
- Income: Ingresos, categorizados en una escala de 1 a 8, 1 = menos de 10,000 USD, 5 = menos de 35,000 USD, 8 = 75,000 USD o mas.
Estas variables ofrecen una visión integral de los factores de riesgo, condiciones de salud y estilos de vida de los participantes, permitiendo realizar análisis detallados sobre la relación entre estos factores y la diabetes.
https://archive.ics.uci.edu/dataset/891/cdc+diabetes+health+indicators
Pregunta 3. Tipo de problema
¿Es un problema de clasificación, regresión o cualquier otro? Explica por qué. Si fuera un problema de clasificación, ¿qué tienes que decir sobre la distribución de los ejemplos de cada clase? Si fuera de regresión, ¿qué tienes que decir de la distribución de valores de la variable dependiente?
El problema planteado en este conjunto de datos se enmarca claramente como un problema de clasificación. Esto se debe a que la variable objetivo, “Diabetes_binary”, es categórica y denota la presencia o ausencia de diabetes en los participantes. En los problemas de clasificación, el objetivo es predecir una etiqueta de clase para cada instancia en el conjunto de datos basándose en las características observadas. La distinción entre tener o no tener diabetes representa una clasificación fundamental entre dos estados de salud distintos, lo cual es un objetivo común en el análisis de datos de salud.
Vamos a ver la distribución de los datos de la variable dependiente “diabetes_binary”.
n<-nrow(diabetesBRFSS2015)
diabetesBRFSS2015 %>%
group_by(Diabetes_binary) %>%
summarize(Frequency = n()/n *100) %>%
ungroup() %>%
ggplot(aes(x = Diabetes_binary, y = Frequency, fill =Diabetes_binary)) +
geom_col(position = position_dodge(), alpha = 0.7) +
scale_fill_viridis_d() +
theme_minimal() +
guides(fill = FALSE) +
ylim(c(0,100)) +
xlab("")+ ylab("Percentage %") + ggtitle("Percentage of Diabetic and Healthy people")
El 50% de nuestros datos son personas sanas y el otro 50% son personas diabéticas. Esto es algo muy importante para tener un modelo robusto con una alta capacidad predictiva.
Vamos a estudiar las variables categóricas con respecto a la variable dependiente para observar los datos y para ver si vemos alguna tendencia.
data_long <- pivot_longer(diabetesBRFSS2015[,c(-5, -15, -16, -17, -20, -21, -22)], cols = -c(Diabetes_binary),
names_to = "Variable", values_to = "Value")
plots <- list()
for(variable in unique(data_long$Variable)) {
data_long$Diabetes_binary <- factor(data_long$Diabetes_binary)
data_filtered <- data_long %>%
filter(Variable == variable) %>%
group_by(Value, Diabetes_binary) %>%
summarize(Frequency = n(), .groups = 'drop') %>%
mutate(Percentage = Frequency / sum(Frequency) * 100) %>%
ungroup()
p <- ggplot(data_filtered, aes(x = Diabetes_binary, y = Percentage, fill = Value)) +
geom_col(position = position_dodge(), alpha = 0.7) +
scale_fill_viridis_d() +
theme_minimal() +
labs(title = variable,
x = "", y = "Percentage") +
theme(legend.position = "none") +
theme(plot.title = element_text(size = 12), axis.title = element_text(size = 10))
plots[[variable]] <- p
}
grid.arrange(plots$HighBP, plots$HighChol, plots$CholCheck, plots$Smoker, plots$Stroke, plots$HeartDiseaseorAttack, ncol=3)
grid.arrange(plots$PhysActivity, plots$Fruits, plots$Veggies, plots$HvyAlcoholConsump, plots$AnyHealthcare, plots$NoDocbcCost, ncol=3)
En estas graficas de barras podemos concluir que casi todas las variables pueden estar muy relacionadas con pacientes con diabetes. Algunas influyen a tener a diabetes como: “HighBP”, “HighChol”, “CholCheck”, “Smoker”, “Stroke”, “HeartDiseaseorAttack”, “DiffWalk” y “Sex” (hombre más propenso). Con esto podemos decir que unos tener una mala salud como una presión arterial alta, colestererol o ser fumador influyen en tener diabetes. El ictus y el ataque cardiaco podría estar relacionado con las variables anteriores por lo que son tambien variables susceptibles de ser predictoras de tener diabetes. En el caso de los chequeo del colesterol en los ultimos 5 años aumenta en pacientes diabeticos. Esto es debido que al ser diagnosticados estos pacientes requieren un mayor seguimiento por parte del personal sanitario.
Hay otros casos que las variables influyen en no tener diabetes como: “PhysHlth”, “Fruits”, “Veggies”. Es decir que hacer ejercicio fisico, comer frutas y verduras ayudan a evitar tener diabetes.
Hay tres variables que no varían mucho entre pacientes sanos y diabéticos que son “HvyAlcoholConsump”, “AnyHealthcare” y “NoDocbcCost”.
Vamos a mirar ahora las variables con varias clases.
# Incluye Diabetes_binary en el pivot_longer, si es necesario
data_long <- pivot_longer(diabetesBRFSS2015[,c(1,15, 20, 21, 22)], cols = -c(Diabetes_binary),
names_to = "Variable", values_to = "Value")
plots <- list()
for(variable in unique(data_long$Variable)) {
# Asegura que Diabetes_binary sea factor si aún no lo es
data_long$Diabetes_binary <- factor(data_long$Diabetes_binary)
# Calcula el porcentaje de Diabetes_binary por cada Value en la variable actual
data_filtered <- data_long %>%
filter(Variable == variable) %>%
group_by(Value, Diabetes_binary) %>%
summarize(Frequency = n(), .groups = 'drop') %>%
mutate(Percentage = Frequency / sum(Frequency) * 100) %>%
ungroup()
# Crea el gráfico
p <- ggplot(data_filtered, aes(x = Value, y = Percentage, fill = Diabetes_binary)) +
geom_col(position = position_dodge(), alpha = 0.7) +
scale_fill_manual(values = c("Healthy" = "#00AFBB", "Diabetic" = "#FC4E07")) +
theme_minimal() +
labs(title = variable,
x = "", y = "Percentage") +
theme(legend.position = "none") +
theme(axis.title = element_text(size = 9))
# Guarda el gráfico en la lista
plots[[variable]] <- p
}
grid.arrange(plots$GenHlth, plots$Age, plots$Education, plots$Income, ncol=2)
Aquí ya vemos algunas distribuciones interesantes. En estos gráficos vemos en color azul los pacientes sanos y en rojo los pacientes diabéticos.
En general, las personas sanas se sienten mejor que los pacientes diabéticos ya que el número 1 en la escala es estado “excelente”. Con respecto a la edad, vemos que la mayoría de personas entrevistadas para este estudio se situan entre los 50 y 72 años (7-11). La diabetes puede clasificarse principalmente en dos tipos: Diabetes de tipo 1 y diabetes de tipo 2, cada una de ellas con diferentes edades típicas de aparición. La diabetes de tipo 1, menos frecuente, suele diagnosticarse en niños, adolescentes o adultos jóvenes, pero puede aparecer a cualquier edad. Por otro lado, la diabetes de tipo 2 es más frecuente en adultos mayores de 45 años, aunque cada vez se diagnostica en personas más jóvenes, incluidos adolescentes y adultos jóvenes, debido al aumento de las tasas de obesidad. Parece lógica esta distribución de los datos ya que el mayor número de muestras se situan en esta franja crítica de aparición de la enfermedad, aunque para saber si la variable edad es un factor de riesgo, se debería tener el mismo número de muestras a cada edad. Considero por tanto que esta variable contiene un sesgo importante para nuestro modelo.
Otro sesgo importante que veo en estos datos, se muestra en las gráficas de distribución de eduacation e income. Vemos que la mayoría de las personas encuestadas para este estudio son personas con alto status socioeconómico ya que tienen altos niveles de estudios y un porcentaje alto de ingresos económicos, sobre todo en la franja de más de 75.000 USD. Al preguntar a este tipo de personas, otras variables se ven muy sesgadas a la hora de optimizar nuestro modelo como por ejemplo: “NoDocbcCost” (no poder costearse un médico), “AnyHealthcare” (tenencia de seguro médico) o “PhysActivity” (actividad física), ya que todas ellas dependen, de manera directa o indirecta, del nivel económico que tenga esa persona.
Vamos a estudiar ahora las variables continuas del dataset.
p1 <- ggplot(data=diabetesBRFSS2015, aes(x=BMI, group=Diabetes_binary, fill=Diabetes_binary)) +
geom_boxplot(alpha=.4) +
scale_fill_viridis_d()+
theme_minimal()+
theme(legend.position = "none",
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
p2 <- ggplot(data=diabetesBRFSS2015, aes(x=MentHlth, group=Diabetes_binary, fill=Diabetes_binary)) +
geom_boxplot(alpha=.4) +
scale_fill_viridis_d()+
theme_minimal()+
xlab("MentHlth (days)") +
theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank())
p3 <- ggplot(data=diabetesBRFSS2015, aes(x=PhysHlth, group=Diabetes_binary, fill=Diabetes_binary)) +
geom_boxplot(alpha=.4) +
scale_fill_viridis_d()+
theme_minimal()+
xlab("PhysHlth (days)")+
theme(legend.position = "none",
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
grid.arrange(p1, p2, p3, ncol = 1)
Pregunta 4. Sesgos
Supón que los posibles sesgos pudieran estar representados por los predictores del dataset. ¿Podríamos determinar visual y numéricamente, con ayuda de PCA si los hay?
En el caso de este dataset hay 19 de 22 variables categóricas por lo que no considero que sea buena idea hacer una PCA para estos datos.
En el apartado anterior comento los sesgos que veo a la hora de recoger los datos y del tipo de personas encuestadas en este estudio.
Pregunta 5. Elección algoritmo
¿Qué algoritmo de entre los que conoces aplicarías según el problema para el que lo has identificado?
Para el dataset diabetesBRFSS2015, que hemos identificado como adecuado para problemas de clasificación, específicamente para determinar la presencia de diabetes, el algoritmo Naive Bayes puede ser una opción atractiva para empezar el modelado.
El algoritmo Naive Bayes es un clasificador probabilístico basado en el teorema de Bayes con la “ingenua” suposición de independencia entre las características. Es especialmente popular en tareas de clasificación de texto pero también puede ser efectivo en datasets de salud como este.
Naive Bayes es conocido por ser un modelo simple y rápido para entrenar lo que lo hace atractivo para datasets grandes o para una exploración inicial de los datos. Funciona bien con variables categóricas lo cual es relevante para el dataset diabetesBRFSS2015, que contiene varias de estas variables, aunque también funciona con variables numéricas.
La suposición de independencia entre las características raramente se cumple en la práctica lo que puede afectar la precisión del modelo. Dado que es un modelo probabilístico puede no ser tan preciso en la clasificación cuando las relaciones entre las características son complejas o altamente correlacionadas.
A menudo considerado como un punto de partida para la clasificación binaria, la regresión logística puede manejar relaciones lineales entre las características y la variable objetivo, pero puede quedarse corta frente a relaciones más complejas que Naive Bayes podría capturar a través de su enfoque probabilístico.
Pregunta 6. Transformación
¿Sería necesario aplicar alguna transformación a los datos? ¿Por qué? Si la respuesta es afirmativa, explica qué tendría de beneficioso hacerlo.
El algoritmo Naive Bayes es menos sensible al escalado de características numéricas que otros algoritmos, por lo que no considero modificar los datos.
Pregunta 7. Algoritmo
Aplica el algoritmo designado en la pregunta 5 y optimiza el estadístico o estadísticos que consideres pertinente para medir su rendimiento mediante la selección de un modelo óptimo. Comenta los resultados
En primer, lugar vamos a generar los datos de entrenamiento y test de los datos de diabetesBRFSS2015 para entrenar el algoritmo de Naive-bayes. Vemos que el porcentaje de paciente sano y diabético en entrenamiento y test se mantienen.
set.seed(1234)
diabetes.TrainIdx.80 <- createDataPartition(diabetesBRFSS2015$Diabetes_binary,
p=0.8, #Genera un 80% para train, 20% para test
list = FALSE, #Dame los resultados en una matriz
times = 1) #Genera solamente una partición 80/20
trainSet.diabetes <- diabetesBRFSS2015[diabetes.TrainIdx.80, ]
testSet.diabetes <- diabetesBRFSS2015[-diabetes.TrainIdx.80, ]
table(trainSet.diabetes$Diabetes_binary)##
## Healthy Diabetic
## 28277 28277##
## Healthy Diabetic
## 7069 7069Dentro del paquete caret encontramos el modelo naive_bayes y los hiperparámetros específicos que se pueden ajustar dependen del tipo de modelo naive_bayes que se está utilizando. Por lo general, naive_bayes es conocido por su simplicidad y tiene relativamente pocos hiperparámetros en comparación con otros modelos de aprendizaje automático más complejos. Sin embargo, hay algunas variaciones del algoritmo que ofrecen diferentes hiperparámetros para el ajuste, especialmente cuando se trata de modelar la distribución de los datos.
usekernel: Este hiperparámetro determina si se deben usar kernels para estimar la densidad de probabilidad de las variables numéricas. Es útil cuando se asume que la distribución de los datos no sigue necesariamente una distribución normal. Los valores habituales: TRUE o FALSE. El valor predeterminado suele ser FALSE, lo que significa que no se usan kernels y se asume una distribución normal (Gaussiana) para las variables numéricas. Si se establece en TRUE, se usan kernels para una aproximación más flexible de la distribución de los datos.
laplace: Este hiperparámetro se utiliza cuando usekernel es TRUE. Se refiere al factor de suavizado (factor de Laplace) aplicado a las estimaciones de densidad del kernel para evitar el problema de la probabilidad cero. Los valores habituales son un valor numérico, generalmente pequeño (por ejemplo, 0, 0.01, 0.1). Un valor más grande aumenta el suavizado.
adjust: También relacionado con el uso de kernels (cuando usekernel es TRUE), este hiperparámetro ajusta el ancho de banda del kernel utilizado en la estimación de la densidad. Los valores habituales son un valor numérico mayor que 0. Un valor de 1 significa no ajustar el ancho de banda, valores menores de 1 lo reducen y valores mayores de 1 lo aumentan. Ajustar el ancho de banda puede ayudar a afinar cómo se modela la distribución de las variables numéricas.
Debido a que computacionalmente naive_bayes no es costoso, he decidido optar por un entrenamiento de validación cruzada repetitiva con 3 repeticiones y 10 particiones. Permito el análisis en paralelo para aumentar el rendimiento.
set.seed(1234)
NBControl.diabetes <- trainControl(method = "repeatedcv",
number = 10 ,
repeats = 3,
returnResamp = "final",
seeds = NULL,
allowParallel = T)En un primer intento de buscar el mejor modelo realizo una parrilla de hiperparámetros con los siguientes valores dando 68 posibilidades combinatorias:
laplace: 0, 0.2, 0.4, 0.6, 0.8, 1
usekernel: T, F
adjust: 0.2, 0.6, 1, 1.4, 1.8, 2.2
set.seed(1234)
y = trainSet.diabetes$Diabetes_binary
x = trainSet.diabetes[,-1]
mygrid.diabetes <- expand.grid(laplace = seq(0, 1, 0.2),
usekernel = c(T, F),
adjust = seq(0.2, 2.2, 0.4))
myfit.diabetes.cv10.rp3.nb <- train(Diabetes_binary ~ ., data = trainSet.diabetes,
method = "naive_bayes",
metric = "Accuracy",
trControl = NBControl.diabetes,
tuneGrid=mygrid.diabetes
)
(myfit.diabetes.cv10.rp3.nb.best <- subset(myfit.diabetes.cv10.rp3.nb$results, laplace == myfit.diabetes.cv10.rp3.nb$bestTune$laplace & usekernel == myfit.diabetes.cv10.rp3.nb$bestTune$usekernel & adjust == myfit.diabetes.cv10.rp3.nb$bestTune$adjust))
El mejor modelo obtenido es laplace = 0, sin utilizar kernel y adjust= 0.2, aunque no haría falta el parámetro adjust ya que no se utilizará kernel.
Como vemos en las siguiente gráficas, casi todas las iteraciones de nuestro mejor modelo de entrenamiento nos da como resultado un accuracy de entorno 0.717 y un Kappa de 0.433. Estos valores son bastante estables en cada una de las iteraciones ya que la distribución de los resultados es muy estable. Los resultados de accuracy, rondan entre 0.7 y 0.73, y los de kappa entre 0.4 y 0.44.
resample <- myfit.diabetes.cv10.rp3.nb$resample
p1 <- ggplot(data = resample, aes(x = Accuracy)) +
geom_density(size = 1.5,color = "#8B1A1A") +
theme_minimal() +
labs(title = "Accuracy",
x = "Accuracy")+
geom_vline(xintercept = median(resample$Accuracy),
linetype = "dashed") +
annotate("text", x = 0.72, y = 5,
label = paste("Accuracy", "\n", round(median(resample$Accuracy), 3)), color = "#8B1A1A", size = 4) +
xlim(c(0.69, 0.73))
p2 <- ggplot(data = resample, aes(x = Kappa)) +
geom_density(size = 1.5, color = "#f1e6b2") +
theme_minimal() +
labs(title = "Kappa",
x = "Kappa")+
geom_vline(xintercept = median(resample$Kappa),
linetype = "dashed") +
annotate("text", x = 0.435, y = 5,
label = paste("Kappa", "\n", round(median(resample$Kappa), 3)), color = "#c0b88e", size = 4) +
xlim(c(0.39, 0.46))
grid.arrange(p1, p2, ncol = 2)
En la siguiente tabla vemos las variables más importantes para nuestro modelo y vemos que la más importante de todas en considerar tener buena salud, seguida de el indice de masa corporal y la tension vascular alta.
import <- varImp(myfit.diabetes.cv10.rp3.nb)
import$importance <- round(import$importance[order(-import$importance[,1]),], 2)
formattable(import$importance, list(
area(col= Healthy:Diabetic) ~ color_bar("lightblue")
))| Healthy | Diabetic | |
|---|---|---|
| GenHlth | 100.00 | 100.00 |
| BMI | 81.88 | 81.88 |
| HighBP | 81.85 | 81.85 |
| Age | 64.85 | 64.85 |
| HighChol | 61.99 | 61.99 |
| Income | 55.67 | 55.67 |
| DiffWalk | 50.38 | 50.38 |
| PhysHlth | 47.41 | 47.41 |
| Education | 38.79 | 38.79 |
| HeartDiseaseorAttack | 31.24 | 31.24 |
| PhysActivity | 30.14 | 30.14 |
| Smoker | 16.41 | 16.41 |
| Veggies | 12.29 | 12.29 |
| Stroke | 11.33 | 11.33 |
| Fruits | 10.04 | 10.04 |
| MentHlth | 9.63 | 9.63 |
| Sex | 7.82 | 7.82 |
| HvyAlcoholConsump | 6.54 | 6.54 |
| CholCheck | 6.05 | 6.05 |
| NoDocbcCost | 2.86 | 2.86 |
| AnyHealthcare | 0.00 | 0.00 |
A continuación, vamos a predecir los resultados de la variable dependiente con el conjunto de test generado anteriormente.
preds.diabetes <- predict(myfit.diabetes.cv10.rp3.nb$finalModel, newdata = testSet.diabetes)
conf.diabetes <- confusionMatrix(preds.diabetes, testSet.diabetes$Diabetes_binary, positive = "Diabetic")
conf.diabetes## Confusion Matrix and Statistics
##
## Reference
## Prediction Healthy Diabetic
## Healthy 5748 4090
## Diabetic 1321 2979
##
## Accuracy : 0.6173
## 95% CI : (0.6092, 0.6253)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.2345
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.4214
## Specificity : 0.8131
## Pos Pred Value : 0.6928
## Neg Pred Value : 0.5843
## Prevalence : 0.5000
## Detection Rate : 0.2107
## Detection Prevalence : 0.3041
## Balanced Accuracy : 0.6173
##
## 'Positive' Class : Diabetic
## La matriz muestra que 5748 individuos fueron correctamente identificados como saludables por el modelo (verdaderos positivos). 4090 individuos saludables fueron incorrectamente clasificados como diabéticos (Falsos negativos). 1321 individuos diabéticos fueron incorrectamente clasificados como saludables (falsos positivos). 2979 individuos fueron correctamente identificados como diabéticos (verdaderos negativos).
Accuracy: El 61.73% de todas las predicciones fueron correctas. Este valor indica la proporción general de predicciones correctas del modelo.
Intervalo de Confianza del 95% para la Precisión: Entre 60.92% y 62.53%, lo que proporciona un rango estimado de precisión del modelo en la población general.
Tasa de No Información (No Information Rate): Es el mayor de las proporciones de las clases de referencia, que en este caso es 50%. Este valor se utiliza para comparar la precisión del modelo; una precisión significativamente mayor que la tasa de no información indica un modelo útil.
P-value de la Precisión Mayor que la Tasa de No Información: Menor que 2.2e-16, lo que sugiere que el modelo es significativamente mejor que una conjetura al azar.
Kappa: De 0.2345, reflejando una concordancia débil a moderada más allá de la coincidencia por casualidad.
Prueba de McNemar Valor-P: Menor que 2.2e-16, lo que indica una diferencia significativa entre los falsos positivos y los falsos negativos.
Sensibilidad: El 42.14% de los individuos diabéticos reales fueron correctamente identificados, lo que sugiere una capacidad moderadamente baja del modelo para detectar la clase “Diabetic”.
Especificidad: El 81.31% de los individuos saludables reales fueron correctamente identificados, indicando una alta capacidad del modelo para detectar la clase “Healthy”.
Valor Predictivo Positivo : El 69.28% de los individuos clasificados como diabéticos eran realmente diabéticos.
Valor Predictivo Negativo: El 58.43% de los individuos clasificados como saludables eran realmente saludables.
Prevalencia de la Clase “Diabetic”: Se asumió un 50%, lo que indica que las clases están balanceadas en el conjunto de datos utilizado para esta matriz de confusión.
Tasa de Detección para “Diabetic”: El 40.66% del total de individuos fueron correctamente identificados como saludables.
Prevalencia de Detección para “Diabetic”: El 30% del total de predicciones fueron para la clase ‘Diabetic’.
Precisión Balanceada: El 61.73%, igual a la precisión general, dado que la prevalencia de las clases es equilibrada.
Aunque el modelo tiene una buena especificidad para detectar individuos saludables, su sensibilidad y valor predictivo negativos son moderados, lo que indica una capacidad limitada para clasificar correctamente a los individuos sanos y para que las predicciones de la clase “Diabetic” sean precisas. El coeficiente Kappa muestra una concordancia débil, lo que sugiere que hay margen de mejora en la capacidad del modelo para clasificar las instancias más allá de lo que se esperaría por azar. La prueba de McNemar destaca una discrepancia significativa entre los falsos positivos y falsos negativos, lo que puede indicar un sesgo en cómo el modelo realiza sus predicciones.
A continuación muestro la matriz de confusión gráficamente para una mejor interpretación.
preds.diabetes <- predict(myfit.diabetes.cv10.rp3.nb$finalModel, newdata = testSet.diabetes)
conf.diabetes <- confusionMatrix(preds.diabetes, testSet.diabetes$Diabetes_binary)
conf.diabetes_tibble <- as_tibble(conf.diabetes$table)
plot_confusion_matrix(conf.diabetes_tibble,
target_col = "Reference",
prediction_col = "Prediction",
counts_col = "n")
Vamos a ver a continuación si este modelo ha incurrido en underfitting.
# Evaluar el rendimiento en el conjunto de entrenamiento
trainPredictions <- predict(myfit.diabetes.cv10.rp3.nb$finalModel, trainSet.diabetes)
conf_train_diabetes <- confusionMatrix(trainPredictions, trainSet.diabetes$Diabetes_binary)
# Evaluar el rendimiento en el conjunto de prueba
testPredictions <- predict(myfit.diabetes.cv10.rp3.nb$finalModel, testSet.diabetes)
conf_test_diabetes <- confusionMatrix(testPredictions, testSet.diabetes$Diabetes_binary)
conf_train_diabetes$overall## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.6122997 0.2245995 0.6082697 0.6163183 0.5000000
## AccuracyPValue McnemarPValue
## 0.0000000 0.0000000## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 6.172726e-01 2.345452e-01 6.092018e-01 6.252951e-01 5.000000e-01
## AccuracyPValue McnemarPValue
## 6.146280e-173 7.107441e-310Al calcular el accuracy de los valores predichos con los datos de entrenamiento y el accuracy de los datos de test, observamos que en los datos de entrenamiento la precisión es ligeramente menor que los datos de test indicando que nuestro modelo ha incurrido en underfitting por lo que deberiamos considerar otro modelo o ver alternativas de preprocesamiento.
En resumen, este modelo predice muchos falsos negativos. En casos medicos como este, un falso positivo o falsos negativos pueden incurrir en graves consecuencias.